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Section  I 


INTRODUCTION 


Use  of  advanced  fiber-reinforced  composite  laminates  has  been  rapidly  growing  in 
structural  engineering,  e.g.,  in  the  design  of  aircraft,  space  vehicles,  automobiles, 
large-span  roof  structures,  etc.  This  is  due  to  the  high  strength/weight  ratio  and  the 
possibility  for  optimal  design  by  tailoring  the  mechanical  properties  of  structural 
components  for  a  specific  application.  Increasing  use  of  composite  materials  in  the 
design  of  high-performance  vehicles  has  attracted  much  attention  to  the  dynamic 
behavior  of  structural  components  under  service  conditions.  Experimental  procedures 
can  provide  information  on  the  real  behavior  of  structures  to  the  designer,  but  cannot 
cover  all  the  design  possibilities.  Therefore,  it  is  important  to  develop  a  general,  as 
well  as  reliable,  analysis  procedure  which  can  predict  the  response  of  composite 
laminates  under  a  variety  of  service  conditions. 

Considerable  research  effort  has  been  devoted  to  the  development  of  analytical 
procedures  for  the  analysis  of  comjxwite  materials.  This  has  resulted  in  a  variety  of 
laminated  plate  theories  and  solution  methods  including,  among  others,  classical  thin 
plate  theory  [e.g.,  Reissner  1961,  Stavsky  1961],  first-order  shear  deformable  theories 
[e.g.,  Yang  1966,  Whitney  1970],  higher-order  theories  [e.g.,  Whitney  1973,1974,  Nelson 
1974,  Lo  1977,  Reddy  1984a]  and  discrete  laminate  theories  [e.g.,  Srinivas  1973,  Sun 
1973,  Pagano  1978,1983,  Seide  1980,  Green  1982]. 
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Classical  thin  plate  theory  (CPT)  based  on  the  Kirchhoff  hypothesis  assumes  that 
the  transverse  shear  deformation  is  negligible.  For  the  analysis  of  laminated  composites, 
it  is  well  known  [Whitney  1969,  Pagano  1969,1970b,  Jones  1970,  Srinivas  1970]  that 
use  of  CPT  leads  to  underprediction  of  the  transverse  deflection,  overprediction  of 
natural  frequencies,  and  higher  buckling  loads.  This  is  attributed  to  the  fact  that  the 
ratio  of  shear  to  Young's  modulus  is  lower  in  most  composite  materials  than  in 
conventional  isotropic  materials.  Also,  the  error  grows  with  an  increase  in  plate 
thick  n“ss. 

This  theoretical  deficiency  of  classical  thin  plate  theory  was  corrected  by  the  shear 
deformable  theory  [Yang  1966]  in  which  transverse  shear  deformation  was  taken  into 
account,  following  Mindlin's  [l95l]  work,  for  the  dynamic  analysis  of  laminated  plates. 
Since  then,  various  shear  deformable  theories  have  been  proposed,  including  higher-order 
theories  in  which  the  power  expansion  for  displacements  contains  terms  of  order  higher 
than  one.  It  has  been  shown  [Whitney  1969,  Srinivas  1970]  that  first  order  shear 
deformable  theory  may  be  adequate  to  predict  global  behavior  of  laminated  plates,  e.g., 
lateral  deflection  or  fundamental  natural  frequency,  but  it  is  not  better  than  CPT  in 
calculating  in-plane  stresses  because  it  does  not  include  the  contributions  of  higher  shear 
mode.s.  Higher-order  theories  lead  to  improved  estimates  of  in-plane  stress  distributions 
and  of  the  flexural  vibration  characteristics. 

However,  the  shear  deformable  laminate  theory,  whether  it  is  the  first  or 
higher-order  theory,  has  two  critical  deficiencies.  The  first  is  its  lack  of  capability  to 
describe  local  deformation  precisely.  Due  to  this,  it  is  difficult  to  avoid  error  in 
calculating  natural  frequencies  as  well  as  in-plane  stresses  around  laminar  interfaces, 
especially,  when  shear  rigidities  of  adjacent  laminae  are  quite  different  [Sun  1973,  Lo 
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1977].  The  other  deficiency  is  the  violation  of  equilibrium  of  the  plate  because  stress 
continuity  at  the  interface  is,  in  general,  not  satisfied.  The  need  to  eliminate  these 
deficiencies  has  motivated  the  development  of  several  discrete  laminated  plate  theories 
[Srinivas  1973,  Sun  1973,  Seide  1980]  in  which  variation  of  anisotropy  in  the 

laminate  is  properly  incorporated.  The  discrete  laminate  theory  not  only  removes  the 

drawbacks  of  shear  deformable  theories  noted  above,  but  it  also  allows  different 

boundary  conditions  to  be  specified  in  each  layer.  It  may  be  regarded  as  the  most 
general  approach  capable  of  accurately  describing  the  mechanical  behavior  of  any  type 
of  laminated  plates.  Use  of  discrete  laminate  theories  appeared  to  give  better  in-plane 

stress  disribution  [Seide  1980]  and  more  accurate  natural  frequencies  [Sun  I97.^j. 
Ho^^ever,  this  theory,  in  general,  involves  a  large  number  of  field  equations,  and 
consequently  makes  the  problems  quite  complicated. 

A  basis  often  used  for  laminate  theories  is  to  assume  a  pattern  of  variation  of 
displacements  over  the  thickness  of  the  plate.  In  such  theories,  which  allow  for  shear 
deformation,  the  constitutive  relations  of  transverse  shear  are,  in  general,  not  satisfied. 
As  a  result,  it  is  not  possible  to  avoid  some  error  in  evaluating  the  laminate  stiffness. 
Since  the  effect  of  transverse  shear  deformation  is  significant  in  laminated  composites, 
accuracy  of  analysis  can  be  considerably  affected.  In  particular,  its  effect  becomes  more 
critical  in  thick  laminates  or  hybrid  laminates  made  of  layers  with  drastically 
different  material  properties.  Many  attempts  have  been  made  to  treat  the  shear 
deformation  realistically,  but  a  standard  procedure  applicable  to  laminates  of  arbitrary 
construction  is  not  available. 

Since  the  boundary  value  problem  of  a  structure  constructed  with  composite 
laminates  is  extremely  complex,  approximate  numerical  techniques  are  often  used  to 
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obtain  the  solution.  The  most  popular  tool  has  been  the  finite  element  method  which 
is  usually  based  on  a  variational  formulation.  Several  different  types  of  element 
geometries,  interpolation  schemes  and  formulation  strategies  have  been  introduced,  (e.g., 
Mawenya  [l974l  Reddy  [1980,1984b],  Bhashyam  [1983],  and  Putcha  [1986]).  To  provide 
the  basis  for  different  possible  formulations,  Al-Ghothani  [1986]  presented 
complementary  variational  formulations  of  the  discrete  laminate  theory  of  dynamics  of 
laminated  plates  following  Sandhu's  [1970,1971,1975,1976]  procedure.  Various  extended 
and  specialized  forms  of  the  general  variational  principle  were  discussed.  However,  he 
failed  to  derive  variational  principles  for  the  direct  formulation  w'hich  provides  another 
and  often  more  useful  approach  for  construction  of  approximate  solution  procedures. 

As  part  of  the  current  research  program,  reliable  procedures  were  to  be  developed 
for  the  analysis  of  stresses  and  deformations  in  delamination  specimens  of  composite 
laminates  allowing  for  the  coupling  of  flexure  and  extension.  This  required 
development  of  a  theoretical  model  which  could  realistically  describe  the  mechanical 
behavior  of  composite  laminates.  The  discrete  laminate  theory  was  selected  as  quite 
general.  This  was  extended  to  include  constitutive  coupling  of  force  resultants  in  the 
lamina.  In  Section  II,  the  field  equations  of  a  discrete  laminate  theory  based  on  the 
assumed-displacement  field  are  summarized  following  Srinivas  [1973],  Sun  [1973]  and 
Seide  [1980],  and  its  somewhat  ad  hoc  treatment  of  transverse  shear  deformation  is 
discussed.  Section  III  presents  a  procedure  based  on  a  generalization  of  Reissner's 
method  to  incorporate  the  effect  of  transverse  shear  deformation  in  a  consistent  manner. 
A  variational  formulation  of  the  consistent  shear  deformable  discrete  laminate  theory  of 
laminated  composite  plates  is  proposed  in  Section  IV.  Direct  as  well  as  complementary 
formulations  are  discussed.  In  Section  V  a  finite  element  discretization  procedure  is 
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introduced.  In  Section  VI,  application  of  the  finite  element  code  to  evaluation  of 
stresses  in  some  cross-ply  and  angle-ply  free-edge  delamination  specimens  is  described 
along  with  an  application  to  free  vibration  analyses. 

The  development  of  the  coupled  shear  theory  discussed  herein  is  an  important  step 
forward  in  obtaining  reliable  estimates  for  stresses  and  deformations  in  laminated 
composites.  Clearly,  the  new  theory  has  certain  limitations  including  its  assumptions  of 
vanishing  transverse  strain.  Further  refinements  on  introducing  coupled  relations  for  the 
other  force  resultants  besides  shearing  forces,  and  allowing  for  variation  of  transverse 
stress  over  the  thickness  of  the  laminate  is  apparently  necessary  for  reliable  estimation 
of  stresses  in  a  composite  laminate. 
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Section  n 


FIELD  EQUATIONS  OF  THE  DISCRETE  LAMINATE  THEORY 
OF  COMPOSITE  PLATES 


2.1  INTRODlJCriON 

In  this  section,  field  equations  of  the  discrete  laminate  theory  for  dynamics  of 
laminated  plates  are  summarized  using  the  kinematic  assumptions  proposed  by  Srinivas 
[1973],  Sun  [1973],  and  Seide  [l980].  The  domain  of  definition  of  all  functions  is  the 
Cartesian  product  Rx[0,oo),  where  R  is  the  closure  of  the  open,  connected  spatial 
region  R  occupied  by  the  plate  and  [0,oo)  is  the  positive  time  interval. 

We  consider  a  laminated  plate  of  uniform  thickness  h  composed  of  an  arbitrary 
number  of  thin  layers,  in  which  each  layer  is  assumed  to  be  homogeneous,  linear 
elastic  with  its  material  axes  not  necessarily  coincident  with  the  geometric  coordinate 
axes  (Figure  1.).  For  the  Cartesiar  reference  frame  used,  the  origin  is  located  in  the 
bottom  surface  of  the  plate  (i,— x,  axes)  with  x,  axis  normal  to  this  plane.  Also,  in 
each  layer  a  local  coordinate  system,  x‘/’,  is  set  up  in  a  similar  way  with  the  range 
of  Xj*’  limited  to  the  thickness  of  layer. 
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(k+1) 


Layer  N 


Layer  : 


Layer  4 


Layer  3 


Layer  2 


Layer  1 


2^  FIELD  EQUATIONS  OF  LINEAR  ELASTOSTATJCS 

Differential  equations  of  equilibrium  for  linear  elastostatics  are: 

where  and  /,  are  components  of  the  symmetric  Cauchy  stress  tensor  and  the  body 
force  vector  respectively.  Here,  and  in  the  sequel,  we  use  standard  indicial  notation. 
Roman  indices  take  on  the  range  of  values  1,  2,  3  and  greek  indices  the  values  1,  2. 
Summation  on  repeated  indices  is  implied  except  where  indicated  otherwise.  The 
superscript  ik)  denotes  the  k'^  layer  and  is  not  summed.  Parenthesis  around  a  single 
index  indicate  ”no  sum”  on  that  index.  Parentheses  around  a  pair  of  indices  denote 
symmetric  part  of  the  tensor  defined  by  the  pair.  Indices  following  a  subscripted 
"comma”  denote  partial  differentiation  with  respect  to  the  spatial  co-ordinate  defined  by 
the  subscript. 

For  small  deformations,  the  strain-displacement  relations  are: 

e  =  —  (u  .  +  u  )  (2) 

>j  2 


where  and  «,  are  components  of  the  strain  tensor  and  the  displacement  vector, 
respectively. 

For  isothermal  elasticity,  the  constitutive  equations  are: 


<7  =  E  ,,e,, 

ij  tjk!  k! 


(3) 


where,  because  of  the  symmetry  of  e,^  and  <r,^,  the  components  E^^,^  of  the  elasticity 
tensor  satisfy  the  symmetry  relation 


E  =  E ....  =  E 

ijk!  jtkl  tjlk 


(4) 


Further,  assuming  the  existence  of  an  energy  function  implies  For  a 

general  anisotropic  uiaterial,  the  elasticity  tensor  with  components  has  21 

independent  elements.  If  inverse  of  (3)  exists. 
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(5) 


e  =  C  .,u.. 

IJ  lfk(  LI 


where  C  are  components  of  the  compliance  tensor. 


23  SPECIALIZATION  TO  A  LAMINATED  PLATE 


2.3.1  Kinematics 


For  a  laminated  plate  subject  to  bending  and  stretching,  in  order  to  reduce  the 
problem  to  one  in  two  dimensions,  the  functional  dependence  of  the  displacements  upon 
the  transverse  coordinate  JCj  is  made  explicit.  Often,  the  in-plane  displacements  are 
assumed  to  vary  linearly  within  each  layer  and  the  thickness  stretch  is  assumed  to  be 
negligible.  Mathematically,  for  the  /c""  layer,  this  can  be  expressed  as 


where  w'*’  are  the  associated  displacements  at  the  bottom  surface  of  the  layer; 
and  are  the  routions  of  the  cross  section  of  the  k'^  layer.  For  small  deformation, 
the  kinematic  relations  for  the  k'^'  layer  are  (2): 

e=—(u+u)=:u,,  (8) 

Substituting  (6)  and  (7)  into  (8),  the  strain-displacement  relations  for  layer  become 

(9) 

=  (10) 

4'  =  o  (11) 


where  the  kinematic  variables  are  defined  as 
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(12) 


-U) 


(*)  _ 


1  /_«)  .  _(X)>.  _  _(X) 
=  “(0.0) 


(13) 


2.3.2  Equilibrium  Equations 

The  three-dimensional  equations  of  motion  of  the  k'^  layer  are 

u )  ,  jn  u)  (t) 

a  +  /  =  p  u 

I ^  r  ^  t 


(14) 


where  p"’  is  the  mass  density.  Here,  superposed  dots  denote  time  derivatives  of  the 
order  denoted  by  the  number  of  dots.  Regarding  a'\f,  /*/’,  uj*’  as  functions  of  x,,  (14) 
is  equivalent  to 


/■(a-'*'  +/‘>-p“V*’);t>,=0 
l  iy.y  ^  t  ^  i  3  3 


for  n=0,  1,  .....  oo.  The  integration  leads  to  a  countable  set  of  equations  involving 
functions  of  x,  and  Xj.  As  an  approximation,  the  values  n  =  0,  1  are  generally  used. 
Fvidently,  higher  order  equilibrium  theory  would  use  higher  order  of  n  as  well.  For 
the  k"'  layer,  integration  of  (14)  and  the  first  moment  of  two  of  the  equations  (viz. 
i  =  l,  2)  over  the  thickness  of  the  layer,  for  the  displacement  assumptions  (6)  and  (7), 
gives  [Al-Ghothani  1986] 

(15) 


^U)  pU)  _  pUtAk)  _  _  Q 

00,0  Or  a  O  0^0 

..U)  Ak)  ^  „(<)  .  .  Ak)  Ak)Ak)  M)Ak)  _  „ 

^00,0  "Qo  +^0  +^^o  “o  *^0 

q'*’  +  (r'^^-r'*  *’)  +  y*’  -  =  0 


(16) 

(17) 


where 
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qV  =  l<l’ 

*o 


dx 


a) 


( =  j  (1. 4*’)  0-1^ 


(  F^'  \  =  )*  (1,  Xj ’)  dx 


a) 

3 


^3  =  /  A  ^^^3 

•?) 

r;*’  =  <r;*>(xf  =  t,)  =  cr!.*^'^(xf ‘^  =  0) 

r‘,-'  =  <r<‘>(x“Uo)  =  <r“->(.<‘-'>  =  ,..,) 

and  tj  is  the  thickness  of  the  *"■  layer. 


U) 

3 


(18) 

(19) 

(20) 

(21) 

(22) 

(23) 

(24) 


2.3.3  Constitutive  Equations 

For  a  composite  lamina  having  material  symmetry  with  respect  to  its  middle 
surface,  coupling  of  the  extensional  stresses  and  the  shear  strains  vanishes  and  (3) 
reduces  to  [Al-Ghothani  1986] 

(25) 


U)  _  ..«)  U)  .  p«)  U) 

"cy0  ^o0y(^y0  ^0033^33 


r‘*J  =  2  £**3  3^  *3 

or3  a3y3  >3 


”33  ^33y6  y6  ^  ^3333  33 


(26) 

(27) 
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Substituting  (25)  into  (19)  and  using  (ll),  the  constitutive  equations  of  bending  and 
stretching  are  obtained  in  terms  of  plate  kinematic  variables  and  force  resultants. 
Namely, 


Ak)  ik) 

yb 

J-U) 

(k) 

«) 

afiyb  ^Ckfiyb 

K  . 
>6 

where 


,  At)  U)  U) 


(«.• 


(2)  3 

J. _  _L) 

2  ’  3 


(29) 


It  is  well  known  from  the  exact  elasticity  solutions  [Pagano  1969,1970a]  that  the 
transverse  shear  stress  distribution  is  close  to  parabolic  over  the  thickness  of  each  layer. 
In  (10),  however,  the  transverse  shear  strains  are  constant  through  the  thickness  of  a 
layer,  which  implies  the  constant  shear  stresses  through  (26).  Furthermore,  if  interface 
continuity  requirement  of  the  transverse  shear  stress  is  enforced,  the  shear  stress 
distribution  becomes  constant  over  the  thickness  of  the  entire  laminate,  which  is  far 
from  the  real  situation.  As  a  result,  direct  use  of  (26)  for  obtaining  the  plate 
constitutive  equations  yields  an  error  in  the  evaluation  of  the  plate  stiffness.  The 
usual  measure  to  avoid  this  error  is  to  multiply  the  shear  stiffness  in  (26)  by  a 
coefficient  K,  but  a  standard  method  to  determine  the  value  of  K  is  not  available. 
Therefore,  a  procedure  to  obtain  the  shear  stiffness  matrix  which  takes  into  account 
parabolic  distribution  of  shear  stress  needs  to  be  developed  so  as  to  ensure  reliability  of 
the  theory.  This  issue  is  addressed  in  Section  III  by  developing  consistent  transverse 
shear  constitutive  relations  which  allow  for  realistic  stress  distributions. 
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2.4  BOUNDAJty  CONDITIONS 


For  the  k*''  layer,  the  boundary  conditions  associated  with  the  field  equations  are: 


on  S**^x[0,oo) 

(30) 

on  S^j^x[0,oo) 

(31) 

on  5**^x[0,oo) 

(32) 

_«)  .A 

u  =  u 

a  ofi 

on  S2^x[0,oo) 

(33) 

on  5*f^x[0,oo) 

4 

(34) 

u)  An,  V 

w  =  tt>  (x^,t) 

on  S*,^^x[0,oo) 

O 

(35) 

where  Xg  are  the  coordinates  along  the  edge  boundary  S**’  of  the  spatial  region  R 
occupied  by  the  plate;  a  circumflex  denotes  the  value  of  the  prescribed  quantity  on 
5^*^:  and  r}g  are  components  of  the  unit  outward  normal  to  S^*\  The  boundary 
segments  5^,*’,  5*/^  are  complementary  subsets  of  S^*\  and  so  are  5^*^  and  Sl*\  Sl*\ 

2J  INlTlAd.  CONDITIONS 

The  initial  conditions  for  the  problem  are 


(36) 

(37) 

w''^.0)  =  wl^\xp 

(38) 

(39) 

(40) 

w<'\x^.o)  =  >v;;^) 

(41) 
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2.6  INTERLAMINAR  CONTINUITY  CONDITIONS 


Since  it  is  assumed  that  all  the  layers  are  perfectly  bonded,  continuity  of 
displacements  and  tractions  along  interlaminar  surfaces  must  be  satisfied.  The 
displacement  continuity  conditions  are: 

(42) 


s'**"  =  s'*'  + 

a  a  K  ^  or 


«•!)  U) 
w  =  w 


and  the  traction  continuity  conditions  are: 


0)/  </)  ,  N  (Ol)/  (X-.1) 

CT  ^  (  = /^  )  =  =0) 


(4.^) 


(44) 


'I’hrough  these  continuity  conditions,  all  the  field  equations  defined  for  each  layer  can 
be  combined  to  give  the  governing  equations  of  the  laminated  plate. 

In  approximate  solution  procedures,  two  distinct  situations  may  arise.  In  case  the 
interlaminar  traction  components  and  the  layerwise  shear  forces  are  admitted  as  field 
variables,  continuity  can  be  directly  enforced.  On  the  other  hand,  if  a  displacement 
type  approach  is  used,  the  shearing  forces  obtained  through  material  contitutive  relations 
can  be  grossly  in  error  if  the  simplistic  kinematic  assumptions  (6)  and  (7)  are  used. 
An  alternative  often  employed  is  to  evaluate  shearing  stresses  from  consideration  of 
equilibrium,  i.e.,  obtaining  through  the  material  constitutive  relations  but 

and  cTj*’  using  (14).  We  discuss  this  point  in  Section  V  where  the  new  theory  is 
applied  to  free-edge  delamination  specimens. 
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Section  HI 


CONSISTENT  TREATMENT  OF  TRANSVERSE  SHEAR 
DEFORMATION 


3.]  INTRODUCTION 

Laminate  theories  based  on  assumed  displacements,  in  general,  do  not  satisfy  the 
constitutive  relations  of  transverse  shear.  Since  the  effect  of  transverse  shear 
deformation  is  significant  in  laminated  composites,  there  could  be  certain  loss  of 
accuracy  in  the  analysis  due  to  this  error.  In  particular,  its  effect  could  be  significant 
in  thick  laminates  and  hybrid  laminates  composed  of  layers  with  drastically  different 
shear  rigidities.  For  this  reason,  to  enhance  the  reliability  of  laminate  theory, 
development  of  a  procedure  to  incorporate  transverse  shear  effect  properly  is  necessary. 

In  this  section,  the  development  of  constitutive  equations  of  transverse  shear  in  a 
consistent  manner  is  described.  The  assumptions  and  notation  of  a  discrete  laminate 
theory  given  in  Section  II  are  used.  The  theoretical  basis  for  development  is  a 
generalization  of  Reissner's  mixed  variational  principle  of  linear  elastic  orthotropic  plates 
to  laminated  composites.  Reissner's  principle  was  stated  on  an  ad  hoc  basis.  Herein,  it  is 
derived  as  an  extension  of  the  general  variational  principle  for  linear  elastostatics  based 
upon  the  general  procedures  for  coupled  linear  problems  introduced  by  Sandhu  and  his 
co-workers  [1970,  1971,  1975,  1976].  A  summary  of  these  procedures  is  given  in 
Appendix  A.  Throughout,  it  is  assumed  that  all  the  functions  are  defined  on 
closure  of  the  open  connected  spatial  region  of  interest  R.  A  rectangular  Cartesian 
coordinate  system  is  used. 
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3.2  COMPLEMENTARY  EORM  OF  FIELD  EQUATIONS  OF  LINEAR 


ELASTOSTATICS 

The  field  equations  (l),  (2)  and  (5)  of  elasticity  can  be  written  as  fol low’s. 


+  <’.13  +  /.  =  0 

(45) 

<7, 

3or.O 

+  <’31.3  +  4  =  » 

(46) 

\-iu  +u  ) 

2  a.P  (3,0 

(47) 

^q3  ~ 

y  ^“o.3  +  “3,0^ 

(48) 

^33  = 

“3.3 

(49) 

C  0  tO"  ,  +  2C  .  ,o-  ,  -f  C 

a^yb  yb  o^y3  >3  o^33  33 

(50) 

C  ,  ,o-  ^  +  2C  ,  ,  or  ,  4-  C 

ot3>o  yo  0f3y3  >3  oi333  33 

(51) 

^33  = 

^33y6^v6  ^^33y3°^>3  ■*"  ^3333^33 

(52) 

Here  we  have  separated  the  equations  involving  spatial  co-ordinate  Xj  from  the  others. 


33  SELF-ADJOINT  FORM  OF  FIELD  EQUATIONS 

Coupled  field  equations  of  linear  elastostntics  (45)-(52)  can  be  written  in 
self-adjoint  matrix  form  as 


0  0 

0  0 


-As  -A 

63  6« 

-L,  0 


0 

As 

03  “> 

2 

A 

A 

0 

03 

dy 

^3333 

2C„  , 

33y3 

^33y6 

2C  ,,, 

4C  ,  , 

2C  , 

o333 

of3y3 

o3y6 

^aP33 

2  Ca(3y3 

^0/3  y6 

U 

y 

-/» 

“3 

-/3 

<^33 

0 

0 

O’  e 

0 

y6 

(53) 


in  which  8.  is  the  identity  tensor. 


ae 


A_ 


and 
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L  =  i(8  +  8„  -^) 

1  2  8« 

L  =  1(8  -14  8,^) 

2  2  ">  a^  ay 

The  operators  on  the  diagonal  of  the  matrix  in  (53)  are  symmetric  tensors.  If  we 
define 

(54) 


</.  8>k  =  f 


dR 


the  off-diagonal  operators  constitute  adjoint  pairs  i.e., 

■^“3’°^33.3^^  ~  ~  '^‘^33’“3.3^  J?  ‘^“3’ ^^33  ^3  ^  0^ 


(55) 

(56) 

(57) 

(58) 


J(  ~  ~  '*'^a3'*^3,a^  X  ^  “s’ ‘^•>3  ^  0^ 

(55)  through  (58)  are  sufficient  to  ensure  self-adjointness  of  the  matrix  of  operators  in 
(53)  in  the  sense  of  (A.25).  Consistent  boundary  conditions  associated  with  the  field 
equations  (53)  are; 


=  u^T)^  on  5, 

(59) 

“3'n3  =  “3^)3  “3")^ 

=  “3’)0  0" 

(60) 

on 

(61) 

-(^33^3  +  0'a3V  =  -  ^ 

on 

(62) 

where  a  superposed  circumflex  denotes  the  prescribed  value  of  the  quantity  over  the 
boundary  surface;  ?,  and  7),  are  the  components  of  the  prescribed  traction  vector  and  of 
outward  unit  vector  normal  to  QR,  respectively.  In  addition,  5,  and  Sj  are 
complementary  subsets  of  QR.  We  note  that  in  a  physical  problem,  each  component  of 
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displacement  or  traction  may  be  specified  over  different  parts  of  the  boundary. 
However,  in  the  interest  of  conciseness,  we  denote  the  part  of  the  boundary  on  which 
displacement  is  specified  as  5,  and  the  portion  on  which  traction  is  specified  as 
This  representation  is  symbolic  and  in  no  way  indicative  of  limitations  on  the  theory 
in  this  respect. 


3.4  A  GENERAL  VARIATIONAL  PRINCIPLE 

Using  the  definition  (A.26),  the  governing  function  for  the  field  equations  (5.'^) 
and  associated  consistent  boundary  conditions  (59)-(62)  can  be  written  as 

fl  ,  =  <u  ,o-  +  <u  ,<T  +  <u„(r  ,  >„  +  2<u  ,/  >  - 

I  o’  o3,3  A  o’  R  3  33,3  R  3  o3.a  R  o  o  R 

+  2  <Uyf  3>^  —  '*'^33’“3,3^Je  ~  ‘^‘^o.3’^“oi.3 ~  “a, 

+  <‘^33-  ^33330-33  +  2033^30-^3  +C33^30-„^>, 

■*'20'a3  ’  ^a333^33  2^or3y3^y3  ^o3y6‘^y6^ Jt 

^0033^33  ^^a0y3^y3  Jl 

+  +  <‘^33*“3"’3-2“3’’3>S, 

+  <tr  ,,  (u  —  2u  )7),>.  +  <(r  ,,  (u.  —  2u.}n  >. 

of3  o  o  '3  S|  of3  ’  3  3  'a  Sj 

-  <u  ,  (o-  „T). +  0-  -7},)  — 2f  >_  —  <u,,  o-,,7),-2t,>^  (63) 

o’  ofi  ^0  o3  '3  0^2  3  33  '3  3  52 

Let  { v}  =  {u^ ,  Uj ,  be  an  admissible  state  corresponding  to  the  set  of  field 

variables  {v}  =  {u„ ,  Mj  ,  o-„g ,  ,  o",,}.  Assuming  that  {v}+XM,  for  \  a  scalar,  is  an 

admissible  state  for  all  X,  i.e.  fl  is  defined  at  every  point  in  a  neighborhood  of  v. 
Gateaux  differential  of  fl,  along  V  is 

A  n  =  — 2<JT  ,,  (u  ,  +  u,  )  —  2C,,  ,cr,-  —  4C  ,  ,o-  ,  —  2C  ,  ,o-  ,>  „ 

t  o3’  0.3  3.0  33o3  33  a3>3  >3  ■y6o3  yl  R 


“"^^33  ’  “33  ~  ^3333^^33  ~  ^^c.333°^o,3  “"  K 

■*‘^‘^“o>  *^03.3  "*"  '*'  '*’  ^33,3  ^o3,o  3^ Jl 

■*‘2'^^a/9’ '*’  ^'^^33' ^“3~“3^3^S, 

+  2  <«y,3.  («3-S3K>S,  +  2  <«^.3.  („-S>3>S. 

-2  < +  <^03^3^  -  ^  >  S,  -  2  <^3  .  (<r 33^)3  +  ‘^.3  V  "  *  3  >  S,  ^^4) 

Because  of  the  self-adjointness  of  the  operator  matrix  (53),  (5l)-(54),  and  linearity  and 
nondegenarcy  of  the  bilinear  mapping,  the  Gateaux  differential  (64)  vanishes  if  and 
only  if  all  the  field  equations  and  boundary  conditions  are  satisfied. 


3-5  EXTENDED  VARIATIONAL  PRINCIPLES  AND  A  SPECIALIZATION 
Equations  (55)  through  (58)  relate  pairs  of  off-diagonal  operators  in  the  operator 
matrix  of  (53)  and  may  be  used  to  eliminate  either  of  elements  in  each  pair  from  the 
governing  function  fl,.  Elimination  of  an  operator  A^J  implies  that  the  state  variable 
need  not  be  in  the  domain  M,j  of  A,j.  This  may  result  in  relaxing  the  requirement  of 
differentiability  of  u^,  thereby  extending  the  space  of  admissible  states. 

Through  this  procedure,  numerous  extended  forms  of  the  function  fl,  are  possible. 
Using  (55)-(58)  simultaneously  to  eliminate  ^^e 

domain  of  the  functional  is  extended  to  include  nondifferentiable  stress  state.  Explicitly, 
this  functional  is 

^2  “  ~2<or^^,  —  2<<r^3,  +  “3,tt^>je  ~2<cr33, 

+  2<u^,f^>g  +  2  <Uyf^>^ 

'^^33’  ^3333*^33 2^33o.3®^«3  "*’^33o/S°^ap^J 
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+  +  2C  ,  ,o-  ,+C  ^,a 

oi3  a333  33  a3y3  >3  or3>6  yo  P 

'*’  ^  '*'  '*'  ^afiyb^yb^  It 

+  2<cr^^.(i/„-ii„)T)^>5^  +  2«T^^,iu^-u^>n^>^^ 

+  2<0-„3.K-^>3>5^  +  2  <0-^3.  (1/3-23)7)^ 

+  2<«^„.t„>,^  +  2<U3.r3>^^  (65) 

This  is  equivalent  to  the  llellinger-Reissner  mixed  variational  principle.  For  this 
functional,  certain  specializations  are  possible  by  constraining  the  admissible  state  to 
satisfy  some  of  the  field  equations.  Assuming  that  (53)^  is  identically  satisfied,  i.e..  the 
constitutive  equation  is  exactly  satisfied  for  the  ’’inplane”  deformations  and  stresses,  ft. 


reduces  to 


ftj  -  ^a.fi^Jt  ■2<<’’o,3'  ^“o,.3  '^'^^33’  “3.3^^ 

+  2  <“.•/•>* +  2  <«,,/,>, 


+  <0’,,.  C„„(r„  +  2C,,  ,cr  ,  +C„  JT  ,>  , 

33  3333  33  33o3  a3  33ofp  op  P 


'^20'„3.  ^a333°’33  2C^33,30'3,3 +C^33^0'3,g>;j 

+  2<o-„p.K-u>p>s  +  2  <tr„.  (u,- 23)113 > 


33’  '**3  “3"''3''S. 


+  2<(T  3,(1/  — 2>q,>_  +  2<cr  (u^  —  u  )-q  >, 
o3  o  o  ^3  .V,  a3  3  3  'o  i 


+  2<u  ,  t  >.  +  2<u,, 

a  a  S2  3  3 


The  only  assumption  to  obtain  ft,  is  that  the  kinematic  and  constitutive  relations  of 
1-2  plane  are  satisfied.  In  connection  with  the  use  of  this  functional  in  deriving  a 
plate  theory  this  point  is  noteworthy  because  most  theories  based  on  the  assumed 


displacement  field  satisfy  this  requirement.  Reissner  [1984]  presented  a  mixed 
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variational  principle  equivalent  to  ftj  'which  was  derived  using  a  Lagrange  multiplier 
technique  and  partial  Legendre  transformation.  For  some  special  types  of  elastic 
materials  with  certain  symmetry  of  material  properties,  the  procedure  for  obtaining  the 
explicit  form  of  the  principle  was  discussed.  However,  an  explicit  expression  of  the 
principle  for  a  general  anisotropic  material  was  not  given.  The  derivation  above  shows 
Reissner’s  ad  hoc  formulation  to  be  a  special  case  of  the  general  variational  principle  of 
linear  elastostatics.  Also,  flj  in  (66)  would  be  more  convenient  than  Reissner's  mixed 
variational  principle  for  the  general  anisotropic  case. 

If  we  assume  further  that  the  displacement  boundary  conditions  on  5,  are 
identically  satisfied,  flj  reduces  to 

fl.=  — <cr,,u  .>  .  —2<(r  (u  -  +i/,  )>.— 2<cr„,u, 

4  or/S  ’  H  of,3  3, Or'  R  33*  3.3  x 

+  2  +  2 

’  ^q333°’33  ^ o.3y  3^  y3  Jt 


+  2<U  ,t  +2<U,,t,>- 
o  ’  a  52  3  ’  3  S2 


(67) 


3.6  A  VARIATIONAL  PRINCIPLE  FOR  A  LAMINATED  PLATE 
The  functional  (67)  written  explicitly  is: 

n.=  |{-^cr.u.  +  a,(M,+u.  )  + 


^03  ^.3  -  y  O’as  ^33  J 


ds 


(68) 


in  w'hich 


e  =  C  ,  ,ar  ,  +  2C  ,  ,  cr  ,  +  C  ,,,  cr„ 

o3  o3>6  >6  ci3y3  >3  o333  33 


(69) 
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6,  =  .cr  ^  +  2C„  ,or  ,  +  (70) 

33  33y6  yb  33>3  y3  3333  33 

Recalling  that  in  the  derivation  of  the  above  functional  the  in-plane  kinematic  and 
constitutive  relation  (53),  was  assumed  to  be  identically  satisfied,  with  some  algebra, 
vanishing  of  the  Gateaux  differential  of  n/u„,tr„a)  along  the  path  (v^,T^gX  gives 

0=  A,  =  /  {  cr  V  -f-  T  ,  +  u^  —2e  ,)  +  t  (u  — e  ) 

^  J  'J  “3  O  3.0,  03  33  3.3  33 


-V  f  ]  dR- 

I  I 


/V, 


Using  (55)-(58)  yields 


0=  A,  .ft  .  =  f  {  — (cr  +f)v+r  (u  .  +  u  —2e  )  +  t  Au  -e  )}  dR 

4  J  **  i  t  o3  a, 3  3.0  o3  33  3.3  33 

-  f  vXt  ~t  )ds  (72) 

i  '  '  ' 

Using  the  notation  defined  in  Figure  1,  the  variational  equation  for  a  laminate 
composed  of  N  layers  is 


.4  *  0 


S,  ,( >  I  0 


U}f  (i)  ,,U)  , 

V  (/  —t  )  dk,  ds 

I  I  3 
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3.7  CX)NS'l  JTUTIVIi  EQUAI’IONS  OF  TRANSVKRSF  SHEAR 


3.7.1  Assumed  Transverse  Shear  Stresses 

In  order  to  use  the  mixed  variational  principle,  developed  in  the  previous 
subsection,  to  set  up  constitutive  equations  for  the  force  resultants,  following  Reisner 
[1984]  we  propose  a  state  of  stresses  in  equilibrium.  The  stresses  are  stated  in  terms  of 
the  force  resultants  as  follows.  Assuming  the  components  to  be  linear  in  Xy  i.e. 


a  „  +  „  j.- 


(74) 


where  and  are  independent  of  x,*’  coordinate,  and  using  the  definitions  of 
force  resultants  it  is  easy  to  show  that 


«)  U) 


(75) 


The  equilibrium  equations  of  elasticity  are,  separating  the  in-plane  equilibrium  equations 
from  the  transverse  ones,  and  ignoring  inertia  terms: 


U  ) 

,  (<) 

+  0-11 
o3.J 

+  , 

Ji) 
^  o 

=  0 

3o,a 

.  U) 

+  ^33.3 

+  , 

M) 

>  3 

=  0 

Integrating  (76)  with  respect  to  Xj*’, 


(76) 

(77) 
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Substituting  (75)  into  (78), 


(*)  J*) 


(*)  Jt) 


^  y,u-i)  _  [4(i^)_3(f3_)2j^a)  _ 

Of  3  Or  ‘  I  ^  ^  ‘  I  *  *  ap,0  or 

k  i  k  A  t 

Substituting  (15)  and  (16)  in  (79)  and  again  ignoring  inertia  terms, 


(/)  U) 

X,  .  .  A-, 


a  '  =  f  +  [4(^)-3(-^)  ][/••  ’  + (7  -r  'O] 

a3  o  ‘  a  o  D 

I  I 

U )  U ) 

f  4  4  O  *  Of  O 

A  A  A 

In  case  of  no  body  force,  i.e.  =  GJ,*‘ =  0,  (80)  reduces  to 


where 


U)  _  ^  M)Ui) 

Of3  *2  Or  ^3  a 


(<•) 

5. 


,1  A,  ,  A, 

C  =  3(^)'  -  4(-^)  +  1 


<3^  =  3(^)^  -  2(^) 

For  a  monoclinic  material,  it  is  only  necessary  to  describe  in  terms  of  t 
evaluate  e„,  in  (69).  Using  the  engineering  notation  for  elastic  constants  Slf  and  5^*’, 
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)  WX )  (X ) 

1  ^55  “45  °^23 

^  -045  244  <^13 


„(/)  „(X)  (i) 

^44  ^45  0^23 

pU)  „U)  U) 
^45  -^55  ^13 


where 


-  O' 

3.7.2  Constitutive  Equations  for  Shear  Resultants 

Neglecting  e,,’  and  noting  that  Uj‘’3  =  0  for  this  formulation,  (73)  yields  the  Euler 
equations  for  transverse  shear 


/ 1  / 


A  X  - 1  0 


Using  (6),  (7),  (81)  and  (82),  and  denoting  the  ’Variation”  in  any  quantity  by  the 
prefixed  symbol  8,  (83)  can  be  rewritten  as 

.V  ^ 

0  =  r  £  { 8(/f^2  ^  +  wj  )  - 

A  t’l  ^5 


where 


[//‘"f  =  [Q^'\  f  ‘I  /’*'*] 


=  IL^’,  4*'.  L^j’] 


(X) 

(X) 

(X) 

1 1 

”l2 

"13 

(X) 

(X) 

n<‘’ 

12 

”22 

"23 

(X) 

(X) 

(X) 

13 

"23 

"33 

25 


«) 


/«)  _  r  M 

■  4  ' 


j=  1.2.3 


ij=  1,2,3 


(88) 


(89) 


nxplicit  evaluation  of  integrals  in  (88)  and  (89)  gives 


1 . 


(/)  6 
577’ 


=  0 


(/>  «)  1 

n,  =  n,  =  —  — 

I’  13  10 


U)  _  (/)  _  2  . 

”22  ”33  15 


u)  _  1  , 

30 


Vanishing  of  the  integral  in  (84),  for  arbitrary  values  of  SQ^  and  S7|7  gives  the 
following  constitutive  equations: 


Jk) 

*^55  ■^45 

6 

1 

rUk-li 

^  1 

1 

rjik) 

*  1 

5^ 

10 

^  2 

10 

rJLk) 

'  2 

*=  1,  2,  .  .  .  N(90) 


J_ 

10 


30 


15 


,«)  ^iOl 

‘ss  ^4S 

ioO)  o<o! 

'•^45  ^44 


„(t)  „U) 
■^55  *^45 

'‘’45  ^44 


^1 


Q 


[yU*-i)i 

li 


10 


^*+1)  ^*+1) 
*^55  ^45 

^*♦1)  ^*+1) 


liLL 

30 


U*+i)  „a*i)] 
•*55  ^45 

1^0+1)  ^a+i)j 

'  45  44 


rrU+l)i 


^i) 

^U) 

5“'” 

0-.) 

55 

45 

+ 

55 

45 

^(O 

^(O 

^U*l) 

UM) 

45 

^44 

45 

44 

r.0) 


*=  1.  2,  .  .  .  N-1 


(91) 


For  a  laminate  of  N  layers,  (90)  and  (91)  constitute  2(27V  — 1)  equations.  These 
equations  may  be  solved  for  7^*’  and  in  terms  of  +  To  do  this,  it  is 
convenient  to  rewrite  these  in  matrix  form  as 
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K«.= 


N. 

0 

t. 

t. 

-— M, 

N, 

_  3 

30  2 

2 

t. 

30 

0 

30  3 

Symm. 


— M  , 
30 


M  , 
30  "  ' 


— 

30  " 


xl  =  [Q\'\di\d'\d'\. . . 

. ,d:\d;) 

(96) 

7  _  r^O  ™(1)  ^2)  ™(2) 

X*  U  ,  .  ^  2  .  4  ,  .  4  2  .  •  •  • 

1 )  »y^(n“  1  ^1 

•  ••••f'l  *  ^  2  ^ 

(97) 

a  ir  j  •  72  ’  41  •  72  ’  •  ’  • 

(n-l)  (n-1)  (n)  ,(o)i 

•  •  •  •»  >1  » y2  ’  ’  ^2  ' 

(98) 

»7  _  1  r  (0)  (0)  „  „ 

~  10  ’  ?2  *  •  •  • 

..0.0,s",g‘”>) 

(99) 

r  _  1  r  (0)  ,  (0)  ^  o 

30  • 

(100) 

M)  M) 

M  = 

^  M)  M) 

*^45  '^44 


Jt=l,  2,  ....  ^  n-l 


0)1  c^-O  oO)  T^*) 

*1  _  ^S5  *45  ^1 

^O)  “  ^t)  ^l)  jit) 


k=  t  or  n 


(t)  ,(t)  . 

V  =  d>  +  w 

»  A  ^  4V  r 


k=  1.  2.  .  . 


(104) 


Solving  (92)  to  eliminate 


28 


(105) 


where 


K  = 

R  = 


aa  ab  bo  ao 


R  -K,KJt^+t 

a  ab  bb  b  a 


Inverting  (105), 

X  =  K  'r  =  AR 

'(i 


(106) 


liquations  (105)  and  (106)  represent  the  relations  between  the  shear  forces  and  the 
shear  strains.  Here  IC  is  symmetric  because  of  symmetry  of  and,  therefore,  A  is 
also  symmetric.  In  (104)  and  (105),  R  depends  upon  the  shear  stresses  specified  on  the 
laminate  surfaces  so  that  the  constitutive  equations  of  transverse  shear  include 
dependence  upon  these  quantities. 

In  general,  K  and  A  are  full  matrices,  '.hus  (105)  and  (106)  may  be  rewritten, 
in  the  absence  of  surface  tractions,  as 


^  =  ZdV’e'.-'’ 


and 


(107) 


(/K 


^=1,  2, 


(108) 


where  and  are  coefficients  defined  by  the  material  properties,  thickness  of 

layers  and  stacking  sequence  of  a  laminate.  From  these  relations,  it  is  seen  that  the 
shear  force  in  a  layer  is  a  linear  combination  of  the  transverse  shear  strains  of  all 
other  layers  and  vice  versa.  This  result  is  due  to  continuity  oi  shear  stresses  in  the 
interfaces  and  shows  that  conventional  approaches  to  handle  transverse  .shear  are  not 
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appropriate.  Also,  symmetry  of  matrices  K  and  A  implies  =  A,'/*'  and 
which  means  that  the  contribution  of  unit  shear  strain  in  the  layer  to  the  shear 
force  of  the  A'*  layer  is  the  same  as  the  shear  force  in  the  layer  caused  by  the 
shear  strain  in  the  layer. 


3.7.3  Specializations  to  the  Mindlin-Type  Laminate  Theory 

The  prcxedure  described  above  can  be  used  to  obtain  the  shear  constitutive 
equations  of  Mindlin-type  plate  theory.  Tor  a  homogeneous  isotropic  plate,  (90)  may  be 
written  as 


G, 

0 

-b 

h 

s- 

T;+r; 

G, 

6 

0 

<^2  +  ^,2 

^2 

T^  +  r- 

(109) 


where  T* ,  (a=l,2)  denote  the  shearing  stresses  specified  on  the  top  and  the  bottom 

surface  respectively.  If  the  plate  surfaces  are  traction-free,  the  relation  (109)  reduces  to 
Reissner's  [1947]  shear  constitutive  equations  with  the  shear  correction  factor  k=5/6. 

For  Mindlin-type  laminate  theory  [Yang  1966,  Whitney  1970],  rotation  of  the  plate 
cross-section  is  constant  and  the  plate  shear  force  resultants  are  the  algebraic  sum  of 


shear  forces  of  all  layers,  i.e.. 


4  - 1 

In  this  case,  the  shear  constitutive  equations  (IO8)  reduce  to 


G  = 


*-i  j-t 


u./y 


(03  +  >V^) 


(no) 


(111) 
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3.7.4  Traction-Free  Edges 


For  the  case  of  free-edge  delamination  specimens,  the  transverse  shear  stress  at 
the  free-edge  is  known  to  be  zero.  This  implies  that  7^*^  and  Q*/’  at  the  free  edge  are 
zero.  Consequently,  the  known  quantities  7^^^  cannot  be  condensed  out  of  (92). 
Explicitly  specifying  7^/’  =  0  and  Q?'  =  0,  (92)  may  be  rewritten  in  the  form 


(112) 


where 


X(i)a  [Qi  *2,  ,  •  •  . 

. . . .  ,e‘;’i 

(113) 

^(2)a  IQ2  '  ^2  ’  •  '  • 

. . . .  ,e';’i 

(114) 

T  _  [•■T’d)  Tr’(2) 

A;  1  i  ,  .  •  •  • 

•  .  .  .  .7’^"''^ 

(115) 

r  _  fy’dl  'T’(2) 

^{2)b  2  ’  2  ’  ’  •  • 

. , 

(116) 

/  1  r  (0)  „  „ 

=  —  « ,  ,0,0,. 

(  1 ).;  10°' 

...  .0, 0,  g'i 

(117) 

I  1  r  .  (0)  rv 

. . .  ..o.o.i^'i 

(118) 

-r  r  (1)  (2) 

=  1^1  ’  Vi  -  ■  •  • 

. yf] 

(119) 

and  and  K,,^  are 

obtained  by  taking  the 

rows  and  columns  corresponding 

to  7^,**  and  Q',*’  from  ,  K,, 

and  respectively. 

Hliminating  from  (112), 


^(l)^(lXi  ^(1) 


(120) 


where 


IT  =  K  —  If  If  *  If^ 

'*■(1)  *^(l)oo  *^(  1)0* 1)06 


If  =Il  —  If  If*T  +T 

“(1)  niXj  '*‘(l)<i6*^(l)W.^l)6^’^<l)o 


In  the  absence  of  surface  stresses,  constitutive  relations  of  the  form  (107),  (108)  at  the 
traction  free  edge  are 


/  i(^)  ,  ii)\  /)  f) 

(0,  +  w, )  =  (2, 

I  1 


and 


(121) 


^=1,  2,  ^  n 


J-i 


where  and  are  the  constitutive  coeeficients  at  the  free-edge. 


(122) 


3.8  AN  EXAMPLE  OF  COUPLED  SHEAR  CONSTITUTIVE  RELATIONS 

For  a  graphite-epoxy  laminate,  made  up  of  12  layers,  each  0.(X)5  inch  thick,  let 
the  material  properties  referred  to  the  material  axes  be 

£,,  =  19.0X10*’  .  £^2  =  £33  =  I.SXIO*’  (pst) 

G,2  =  G,3  =  0.8  X  10*  ,  G^j  =  0.528  X  10*  (psi)  (123) 

^X2  =  *'13  =  •  *'23  =  0.42 

To  study  the  role  of  coupling  in  constitutive  relations  for  shear  forces,  we  consider  the 
stackings  [Oj/903],  and  [ +453 / -45, ],. 
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Table  1  shows  the  transverse  shear  stiffness  coefficients  for  the  [O3/9O3]  laminate 
and  Table  2  contains  those  for  the  [+4537-453]  laminate.  Only  the  coefficients 
corresponding  to  the  transverse  shear  stress  resultants  have  been  listed,  i.e.,  the 
For  the  first  laminate,  Qi*'  and  are  uncoupled;  i.e.,  =  0.  For  the 

second  laminate  the  coefficients  for  and  Qjik)  are  identical  due  to  the  fibre 

orientation  of  45°;  i.e.,  X*,*''’ =  X^*^’.  The  diagonal  terms,  X*a,  represent  the  shearing  force 
in  each  layer  due  to  unit  shear  deformation  of  the  same  layer.  The  off-diagonal  terms 
represent  the  coupling  between  layers.  As  is  evident,  for  the  cases  studied  the 
interlayer  coupling  is  not  "strong”  i.e.,  the  shearing  force  in  any  layer  is  not 
significantly  influenced  by  the  deformation  of  the  others.  Also,  the  effect  is  localized 
i.e.,  the  contribution  of  deformation  of  any  layer  to  the  shearing  force  in  another 
decreases  sharply  with  the  distance  between  the  layers.  Table  3  and  Table  4  show  the 
inverse  of  the  stiffness  coefficients  i.e,  the  compliance  coefficients 

It  should  be  noted  that  in  the  case  where  G,j  =  G23  there  will  be  no  coupling 
between  Q',*’  and  Q'/’  by  virtue  of  being  zero.  Moreover,  the  inter-layer 

coupling  will  be  independent  of  the  fibre  orientation  as  and  ^  will  no  longer  be 
affected  by  the  orientation  of  the  fibres. 


33 


Table  1 


Transverse  Shear  Stiffnesses  for  [O3/9O3]  Laminate 


Layer 

Stiffness  Coefficients 

(k) 

^11 

x‘‘^ 

^11 

x“ 

^11 

^11 

^11 

1 

3447 . 698 

133.869 

22 . 184 

3.137 

0.538 

0.092 

2 

133  869 

3603 . 755 

155.289 

21.958 

3.767 

0  646 

3 

22  184 

155.289 

3576.224 

128,611 

22.066 

3.786 

4 

3.137 

21.958 

128.611 

2404 . 119 

110.513 

18.961 

5 

0  538 

3.767 

22.066 

110.513 

2382.900 

106 . 872 

6 

0.092 

0.646 

3.786 

18.961 

106 . 872 

2382  300 

7 

0.016 

0.111 

0.649 

3.253 

18.337 

106  770 

8 

0  003 

0.019 

0.111 

0.559 

3.149 

18.337 

9 

* 

0.003 

0  019 

0.099 

0.559 

3.253 

10 

t 

0.004 

0.019 

0.111 

0.649 

11 

» 

3(- 

0.003 

0.019 

0.111 

12 

’ 

- 

* 

0.003 

0.015 

*  denotes  coefficients  smaller  than  10 


in  magnitude . 
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=  0. 
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Table  3 


Transverse  Shear  Compliances  A  for  [  O3  /  9O3  ],  Laminate 


Layer . 

Compliance 

Coefficients  ( x  10 

k.2 

k.^ 

14 

k.^ 

kt 

(k) 

M,, 

Mu 

Mu 

Mu 

Mu 

Mu 

1 

0 . 2904 

-0.0107 

-0.0013 

-0.0002 

* 

2 

-0.0107 

0.2784 

-0.0119 

-0.0018 

-0.0002 

3 

-0.0013 

-0.0119 

0  2807 

-0.0148 

-0.0018 

-0.0002 

4 

-0.0002 

-0.0018 

-0.0148 

0.4176 

-0.0191 

-0.0024 

5 

♦ 

-0  0002 

-0.0028 

-0.0191 

0.4214 

-0.0186 

6 

* 

%■ 

-0  0002 

-0.0024 

-0  0186 

0.4214 

7 

A. 

-0.0003 

-0.0023 

-0.0186 

8 

> 

»■ 

ic 

-0.0003 

-0.0024 

9 

« 

* 

-0.0003 

10 

» 

11 

4 

* 

> 

12 

• 

» 

' 

» 

'  denotes  coefficients  smaller  than  10  **  in  magnitude. 


Transverse  Shear  Compliances  Mii  lor  l-^453/-45,  ]  Laminate 
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dMOtas  ce«ff icirata  aaallar  than  10  *  in  aagnituda. 


3.9  DISCUSSION 


For  determining  the  constitutive  equations  for  transverse  shear  in  discrete  laminated 
plate  theory,  a  mixed  variational  principle  of  linear  elastostatics  has  been  derived.  The 
basis  for  derivation  was  the  method  proposed  by  Sandhu  [1970,1971,1975]  for  the 
variational  formulation  of  linear  coupled  problems  with  multiple  field  variables.  The 
variational  principle  is  equivalent  to  Reissner's  mixed  variational  principle  [19841  but 
more  convenient  for  application  to  a  material  with  general  anisotropy.  Using  this  mixed 
variational  principle,  a  procedure  to  obtain  the  constitutive  relations  for  transverse  shear 
has  lieen  developed  for  a  discrete  laminate  theory  which  is  based  on  the  assumptions  of 
linear  in-plane  displacements  and  parabolic  transverse  shear  stresses  over  the  thickness 
of  each  layer.  The  procedure  allows  for  the  interlaminar  continuities  of  stresses  and 
displacements.  Resulting  constitutive  equations  show  that  the  shear  force  resultants  of  a 
layer  are  coupled  with  the  shear  strains  of  the  other  layers  as  well  as  of  different 

directions  (x,  and  Xj).  As  indicated  by  earlier  investigators,  the  shear  stiffness  of 

X,— Xj  and  Xj— x,  sections,  in  general,  are  different  and  vary  with  stacking  sequence  of 
a  laminate.  Also,  the  consistent  shear  constitutive  relations  for  the  Mindlin-type 

laminate  theory  have  been  derived  as  a  special  case.  Actual  computation  of  the  shear 
stiffness  requires  inversion  of  a  square  matrix  with  constant  elements.  This  can  be 
carried  out  with  a  high-speed  digital  computer  without  much  difficulty.  The 

constitutive  relation  for  shear  at  the  freee-edge  or  surfaces  on  which  shearing  stresses 
are  specified  is  different  from  the  relation  when  these  stresses  are  not  specified.  The 
example  of  a  12-layer  graphite-epoxy  laminate  was  considered  using  [0j/90j],  and 
[+45,/— 45j,  stackings  and  the  extent  of  coupling  studied. 
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Section  IV 


VARIATIONAL  FORMULATION  OF  DISCRETE  LAMINATE 

THEORY 


4.1  INTHODUCnON 

Procedures  for  obtaining  approximate  numerical  solutions  to  boundary  value 
problems  are  often  based  on  variational  formulations.  Tor  systematic  development  of 
variational  principles  governing  linear  and  certain  nonlinear  problems,  general  procedures 
have  been  developed.  Mikhlin  [1965]  set  up  the  problem  in  an  inner  product  space 
and  stated  the  basic  variational  theorem  for  a  self-adjoint  boundary  value  problem 
with  homogeneous  boundary  conditions.  For  deriving  variational  principles  governing 
initial  value  problems,  Gurtin  [1963,1964]  used  convolution  product  as  the  nondegenerate 
bilinear  mapping  and  explicitly  included  non  homogeneous  initial  and  boundary 
conditions  in  the  formulation.  Sandhu  [1970,1971]  extended  these  ideas  to  the  general 
linear  coupled  problem.  In  the  context  of  application  of  the  finite  element  method, 
Prager  [l968]  included  in  the  variational  formulation  jump  discontinuities  which  may 
exist  across  interelement  boundaries.  By  introducing  the  concept  of  boundary  operators 
consistent  with  the  field  operators,  Sandhu  [1975]  examined  the  general  case  of  linear 
operators  with  nonhomogeneous  boundary  conditions  and  internal  jump  discontinuities. 

For  mechanics  of  the  fiber-reinforced  composite  laminates,  little  work  has  been 
done  on  variational  formulation  of  the  problem.  Al-Ghothani  [19861  following  Sandhu 
[1970, 1971, 1975,19761  presented  a  variational  formulation  of  dynamics  of  laminated 
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composite  plate.  General  variational  principle  was  derived  based  on  the  complementary 
form  of  an  extension  of  Seide’s  [1980]  discrete  laminate  plate  theory  to  include  inertial 
force,  allowing  for  non  homogeneous  boundary  conditions  and  internal  jump 
discontinuities.  Various  extended  and  specialized  forms  of  the  general  variational 
principles  were  discussed.  However,  he  failed  to  derive  direct  variational  formulation 
which  gives  other  types  of  variational  principles.  Furthermore,  the  laminate  theory  used 
did  not  treat  the  effect  of  transverse  shear  deformation  adequately.  This  effect  is 
important  in  studying  IcKal  deformation  and  pos.sibly  in  modelling  higher  vibration 
mcxJes. 

In  this  section,  a  variational  formulation  of  the  problem  of  vibration  of  a 
laminated  composite  plate  allowing  for  nonhomogeneous  boundary  conditions  as  well  as 
internal  discontinuities  is  developed.  The  theory  used  is  the  one  described  in  Section  II 
along  with  the  constitutive  equations  of  the  transverse  shear  derived  in  Section  III. 
Even  though  the  principal  concern  of  the  research  program  is  the  behavior  of  free-edge 
delamination  specimens  under  static  loading,  inertia  effects  and  arbitrary  geometric 
configuration  and  loading  are  included  in  the  variational  formulation  because  of  the 
ease  with  which  such  generality  could  be  introduced.  Both  the  direct  and  the 
complementary  forms  of  the  field  equations  are  considered.  Extended  variational 
principles  based  on  self-adjointness  of  the  operator  matrix  are  introduced  along  with 
several  specializations.  One  of  their  specializations  formed  the  basis  of  the  finite 
element  approximation  described  in  Section  V. 
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4.2  INTEGRAL  FORM  OF  IT  ELI)  EQUATIONS 

4.2.1  General  Comments 

To  set  up  the  function  governing  the  motion  of  laminated  plates,  it  is  necessary  to 
write  the  field  equations  in  a  way  that  the  operator  matrix  is  self-adjoint  in  a  certain 
space.  The  self-adjointness  of  operators  is  not  an  absolute  notion,  but  rather,  it  is 
relative  to  the  choice  of  an  appropriate  bilinear  mapping.  Thus,  there  are  two  possible 
ways  to  set  up  variational  principles  governing  the  problem;  one  is  to  find  a  bilinear 
mapping  that  makes  the  field  operators  self-adjoint,  and  the  other  is  to  transform  the 
field  equations  so  that  they  can  be  self-adjoint  with  respect  to  a  familiar  form  of 
bilinear  mapping.  For  various  intial-boundary  value  problem,  Gurtin's  [1963,1964] 
procedure,  which  belongs  to  the  latter  approach,  has  been  successfully  applied  [Sandhu 
1971,1987],  and  we  follow  it  for  the  present  problem  although  other  forms  can  be 
used.  Transformation  of  the  differential  form  of  the  field  equations  to  the  equivalent 
integral  form  is  done  by  applying  Laplace  Transform  and  taking  inverse  after 
appropriate  rearrangement.  The  procedure  removes  the  time  derivatives  from  the 
equilibrium  equations  and  includes  initial  conditions  explicitly.  For  the  field  equations 
given  in  Section  II  along  with  the  constitutive  equations  of  transverse  shear  derived  in 
Section  Ill.  integral  form  of  the  field  equations  is  presented  below.  Throughout,  an 
asterisk  (*)  denotes  the  convolution  integral,  i.e. 

r 

u*v=  r  u(t) v(t  —  t) dr  (124) 

We  note  that  the  convolution  satisfies  distributive,  associative  and  commutative  laws. 
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4.2.2  Kinematics 


Equations  (9)  and  (10)  upon  taking  convolution  with  the  time  variable  become: 


ap  ap  3  ap 


ot3  9  ,a 


where,  by  (12)  and  (13), 

+u^^) 

nft  2  "•f 


i  X  n  -  ^  ^  (0  «  +0«  > 

(^P  y  ~(^,p  ~  p,o 


(125) 

(126) 


(127) 

(128) 


4.2.3  Equilibrium  Equations 

Equations  (15)  through  (17)  upon  convolution  with  t  and  appropriate  integration  to 
eliminate  derivatives  with  respect  to  time  give; 


t*N 


t*M 


.it) 


a  ^  ■ 
k  o 


where 


+  i*(7^*'  -  +  z'^'^  =  0 

*'0.0  3  J  3 


O  OO  OO  Oo  oo 


O  OO  OO  OO  OO 


=  f‘^^(wJJ’  +  i*>Vq') 


Or  ~0  O 

(129) 

jf(*Ui)  _  ^  yr«)  _  Q 

a  ^Of  Or 

(130) 

pM^U  )  ^  ^U  )  _  Q 

(131) 

(132) 

(133) 

(134) 


The  initial  conditions  (36)  through  (41)  appear  explicitly  in  the  equilibrium  equations 
above. 
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4.2.4  Constitutive  Equations 


Equation  (28)  upon  convolution  with 


t  yjelds; 


•» 

II 

Ak) 

(*) 

yb 

„{k) 

ik) 

^a^yb 

y6 

and 


(108)  upon  convolution  with  /,  noting  (10)  gives 


’q'"  =  21*  £  x'';’*';; 


The  inverse  relations  are 


_(/) 

-Al. ) 

— ik ) 

yb 

^0,0 

=  /* 

^0/5  y& 

^in^yb 

0) 

■j^k) 

^afiyb 

— ii ) 
^afiyb 

yb 

and 


(135) 


(136) 


(137) 


(138) 


4.2.S  Boundary  Conditions 

As  with  the  field  equations,  the  boundary  conditions  (30)  through  (35)  upon 
convolution  with  t  give 


1*  V 'ti.  =  t*F} 

ofi  ’(S  o 

M) 

on  Sj 

(139) 

M) 

on 

(140) 

^’Or  **0 

on 

(141) 

t*u  =  t*u 

O  O 

on  5^*^ 

(142) 

“a  ~a 

on  5^;’ 

(143) 

t*W  =  t*W 

M) 

on 

(144) 
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4.2.6  Interlaminar  Continuity  cf  Displacements 


Equations  (42)  and  (43)  upon  convolution  with  l  give 


t*u  =t*u  +t*r,d> 

a  Ot 


(145) 

(146) 


43  DIRECT  VARIATIONAL  FORMULATION 


4.3.1  Self-Adjoint  Form  of  the  Field  Equations 

The  field  equations  of  fiber-reinforced  composite  laminated  plate  expressed  in 
integral  form,  (125)-(136)  and  (145)-(146)  can  be  written  in  the  self-adjoint  matrix 
form  as 


A.  B.  D,  ,  0  D. ,  0  •  •  •  •  D,  .  0  D. 

U, 

1  1  1.2  1,3  l,n-l  l,n 

*^1  ~0 

O 

O 

o 

o 

o 

o 

o 

Si 

0 

A,  B,  D,,  0  •  •  •  •  D,  ,  0  D- 

2  2  2.3  2y>- 1  2/1 

^2 

-*•2 

O 

o 

e 

o 

e 

Sz 

0 

A-  B  •  •  •  •  D,  ,  0  D, 

3  3  3m- 1  3m 

^3 

-*•3 

0  0  0  0 

~3 

0 

A  ,  B  ,  D  , 

ri  1  n-  1  n- 

U  1 

1 

0 

"h-  I 

0 

A 

" 

U 

n 

n 

(147) 


or,  symbolically, 

M{y}  =  {Z} 


Here,  only  the  operators  in  the  upper  triangular  region  have  been  entered.  The  below 
diagonal  operators  are  adjoints  of  the  above  diagonal  operators,  i.e.,  the  operator  A,^  is 
adjoint  of  in  the  sense  of  the  bilinear  mapping  used  to  set  up  the  variational 
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formulation.  To  satisfy  the  self-adjointness  condition  (A-25)  in  Appendix  A,  it  is 
sufficient  that  the  elements  of  the  matrix  X  be  adjoints  of  elements  X for  i^j 


and  the  diagonal  elements  be  self-adjoint.  Explicitly,  the  symbolic  operators  applying 
in  (147)  are; 


= 


a 

0 

0 

0 

0 

0 

0 


r*r, 

1  ay 


C*) 


a0y6 


fB] 

0 

0 

0 

0 


Hyfi 


0 

0 

0 

0 

0 

0 

0 


0 

0 

/“’fi 

u 

0 

0 

0 

-/• 


-/*Ba8y6 


0 


0 

0 

0 

/T 


t*D 


,(*) 


aQyb  ^ 

t*  0 


0 

0 

0 

0 

0 

0 

0 


0 

“ray 


0 

0 

0 

-t* 

0 

0 

rs 

orgy 

f 


(148) 


bI  = 


\t*  0  0  t*t^  0  0  0  0  0] 
000  0  00  t*0  0| 


(149) 


<  = 


-1*  0  0  0  0  0  0  0  0 
0  00000  -r*  00 


(150) 


D  = 

‘■j 


|0  000000  0  0 
0  0  0  0  0  0  0  0  0 
0  0  0  0  0  0  0  0  01 
0  0  0  0  0  0  0  0  0| 
0000000  0  0 
0000000  0  0 
0000000  0  0| 

0  0  0  0  0  0  0  -t*  X'^  oi 

ya 

jo  000000  0  01 


ij=  1,  2, 


(151) 


.,T  r_(*)  _(X)  .,(i)  At)  «)  ,,a)  U)  _  «■) 


(152) 
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(153) 


j_r 


=  [7 


,(o 


^U)] 


*  7  /  ”y  y 

-r  ^  0.0,0.  0,  0,  t*Tf,  0.  0] 

-r  ^  0, 0,  ,  0, 0] 


r, 


-^) 


r,  =  i(8  -4 
2  ">  as 


+  s 


o^ 


(154) 

(155) 

(156) 

(157) 

(158) 


Here,  7’|®’  and  7’]"'  are  specified  shear  stress  components  on  the  top  and  bottom  surface 
of  the  plate  and  S^^  is  Kronecker's  delta.  Operator  matrix  in  (147)  is  self-adjoint  in 
the  sense  of  (A-25)  if  the  bilinear  mapping  is  defined  as 

<u,v>  =  r  u*vdv  (159) 

which  is  linear  and  nondegenerate.  Nondegeneracy  of  this  bilinear  mapping  was  shown 
by  Curtin  [l963,1964j. 

The  boundary  conditions  consistent  with  the  operator  equations  (147)  are 


t*  V* ’t).  =  t*fl 

Or^  '0  O 

on  S,  (x  ) 

1  O 

(160) 

t*  =  t*  ^ 

op  'p  o 

on  5,  (X  ) 

3  0. 

(161) 

t*Q,  \  =  t*Q^  X 

eU)/ 

on  S,  (x  ) 

5  Or 

(162) 

t*u^\  = 

on  ^Xz^sufikH) 

(163) 

on  S''W*') 

4  o 

(164) 

71  =  t*  W  71 

'o  'a 

on  (x  ) 

6  a 

(165) 

and  the  internal  jump  discontinuity  conditions  are 
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on 

It  Ck 


on 

3l  Or 


,.(Q«  V  =  ‘’igsXx 


on  s'Ax'*') 

5(  0) 


on 

2t  O 


o«V  (''K 
on  S.  (x  ) 

4/  o 


)'  =  t*(g^^)r) 


Mh  U\ 

on  (x  ) 

6/  a 


The  following  relations  are  satisfied  by  off-diagonal  operators. 

■''“a  ^  jf(k) 

^  a  ’*  ^^ocBr'B'^^k)  ^  ^^^0,0'  *  “o.  ^/9  jCO 


-  <sl*’.<*<w'>P'>j,)+  <<*;.<‘<a“V>a., 

*2f 


*(^) 


<0o  .  ^  ^00^0^  ^k)  ■*■  ^^a0'  ^ 


-  <<‘’-  +  “=<r  '•<*“  V>s'*' 

•^3/  *^4/ 


»(4:) 


<w  ,  r*0  >  ,,  =  —  <f*  w  ,  Q  >  /o 

^0,0  p'>^>  .a  ’  nU) 


.  U)  ,*^U)  .  .  **  (^)  ^ 

-<>V  +  <C„  .f*>v 


-  <>v^*\  +  <cr. 
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4.3.2  A  General  Variational  Principle 

Using  (A-29),  the  governing  function  for  the  operator  equations  (147)  is  defined  as 

n  n- 1  n 

n  =  £  <u,'',  A,u,>  +  £  <  u[.  >  ,„  +  £  <  uj,  cs,_,  > 


X-l  j=l  x-i 


+  L  <=[.  b;u,  +  £  <<•  C"u„,  > 


+  2Z  < U," ■  f,  >  -  2 <  +  2  <  < .  2.  > y., 

/  -  I 

•^Boundary  Terms  +  Internal  Jump  Terms  (175) 

Substituting  (l48Ml7l)  into  (175),  the  explicit  form  of  the  governing  function  is 
obtained. 


'*■'^^0^’'^  •^C»p>6«y6  ^  ^ afi  *  ^ofiyff^yb  ^ /k) 


+  +t*e^'‘  > 


+  <0  . -R  u  —I<p  — 

O  *^0  «P-f3  Of 


+  <2e„,.-2f»X„^ep3  +  f»(2„  >  ,„ 
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+  <Q  ,-i*<P  — «*w  +2/*^  ,>  ,,, 

■^o  ,o  o3 


+  r  {  <u*\  r*r'*^>  +  <w^*’,  r*7’'*'>  } 

Of  ’  a  j^UJ  ^Of  ’  k  a  ’  3 


u) 


o  o  k^ o  jkk)  3  *  piX)  * 


•«) 


+  £(  <?!'’. -»‘2"''’>.„, +  <7^’. ->•»>“•'’>  ,„) 


+  22;^{  <M«  >7*,+  <0^  .i*G^  +y^  >  ,,)+  <w  ,t*F^  +z  >  ,} 


U) 


A±  T?(^^  I 


-2<u‘\  (,)- 2<w‘‘\f*7^f> 

a  o  jj(U  ’  3 

+  2<u''\  f*7^"’>  ,  ,  +  2<«^‘"\  t*t  ,  +  2<w^"\  t*T^f> 

Or  O  ^Or  no  gkft)  3  i 


+  <w-,  r*(er-2CH  +  <<;.  tHu:>-2U-\ 


+  <w"’,  ^*((q1'\,)'-2(/*\7)J  >  +  t*((u“\)'-2(g‘,^\T]p)  > 
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/<)  **//..<<)  V  -»/ _<*)\ 


This  function  is  defined  over  the  set 

f(k\_f  U)  ,(k)  ,(*)  .,u)  „  ^a)  M)  (k)  (k)  ,jXk)  ™a)i 

where  each  of  the  functions  in  the  set  is  sufficiently  smooth  for  the  governing 
function  to  exist.  This  requires  u,!^\  <i>^\  to  be  continuous  and  such 


that  their  derivatives  have  a  finite  number  of  discontinuities.  The  collection  of  all 
possible  sets  within  the  domain  of  fl  is  the  set  of  admissible  states,  l^t 
i-«)i  ,  Ak)  ^U)  -U)  Ak)  .  Ak)  Ak)  Ak)  Ak)  Ak)  Ak)-, 

be  an  admissible  state.  Assume  v  +  Xv€  the  set  of  admissible  states  for  all  values  of 
the  scalar  A,  i.e.  the  function  fl  is  defined  in  a  neighborhood  of  v.  The  Gateaux 
differential  (A-IO)  of  the  function  fl  defined  by  (176)  along  the  path  using  (172) 


through  (174)  to  eliminate  ^1* la* ’^!*ap/A<. •  ^a!s* ^1* !»  “ 
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<-l 


4*> 


+  <mj^.  2«»((</)^j’T,0)'-(g‘^'\T)p>^,,.,  +  2«»((w'^’7)J'-(gJJ>^„ 


M) 


UK 


Ik) 

2/ 


(177) 

The  Gateaux  differential  vanishes  if  and  only  if  all  the  field  equations  along  with  the 
boundary  conditions  and  the  internal  discontinuity  conditions  are  satisfied  because  of 
linearity  and  nondegeneracy  of  the  bilinear  mapping. 

Rearranging  the  terms,  the  governing  function  (176)  may  be  written  as 
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n  1 


a  Of  Ar^or  o  3 

i  =  I 

+2j'{  <u*\  +  ,.,  +  <0‘*\  +  +  <w‘*\t*y*^  +  z‘*^>  } 

o  o  a  ^or  *  ot  o  ’  3  jl^*) 


X-l 


-2<a‘'\  .  -  2<w‘'\ 

o’  o  p(  1 )  ’  3 


>(1) 


+  2<u"\  t*T^"^>  ,  ,  +  2«t>"\  t*t  7’‘"'>  ,  ,  +  2<w‘"\  t*T^"'’>  .  . 

o  ’  o  pyn)  ♦  «  o  j^Uif  ’  3 


«) 


,u)  «) 


+  <  W  ,t*(Q^  7)^-  2Q^  T)J  >  .  t*  (u^  T)^  -  2u^  T)^)  > 


.(*)_ 


!v(i)_ 


t-l 


£{  <<*>,  <*((W>/-2(,“>).)>^..  +  <€■  '•««’)■ -2(J,‘V>^.. 

■^3* 

+  <w<‘',l>((e"r,.)--X,"').„.)>  +  <w‘“,<-((s“V-2(s“Hf>^.l 

■^s#  *^21 

+  <«'".  >,.«  +  <€■  '•«“’"\)’-«sl‘’)».)  > 

(178) 


^*) 

■’*( 


4.3.3  Extended  Variational  Principles 

Equations  (172K174)  relate  pairs  of  operators  in  the  operator  matrix  of  (147). 
These  relations  may  be  used  to  eliminate  either  of  f^ao^  or  2*^*^  either  of 
or  and  either  of  or  w^’  from  the  governing  function  ft  in  (176), 

leading  to  numerous  different  forms  of  variational  formulations.  Elimination  of  an 
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operator  implies  that  state  of  variable  need  not  be  in  the  domain  of 
Where  A,^  are  differential  operators,  this  results  in  relaxing  the  requirement  of 
differentiability  in  u^,  thereby  extending  the  space  of  admissible  states.  In  the  context 
of  the  finite  element  method,  it  is  clear  that  the  extension  of  the  admissible  space 
provides  greater  freedom  in  selection  of  approximation  functions.  To  llustrate  the 
procedure,  we  present  six  extended  variational  principles. 

Using  (172)  to  eliminate  (178)  gives 
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+  -2  <«“>.  >  +  «A;>.  > 
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(179) 

where  need  not  be  differentiable.  In  addition,  eliminating  M^^h,  from  (179)  O, 
reduces  to 
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Also, 
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Q‘^*’  can  be  eliminated  using  (174)  from  (180)  to  give 
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where  none  of  the  stress  resultants  M^Jg  and  d,*’  need  be  differentiable. 
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Alternative!};,  extended  formulations  which  do  not  have  derivatives  of  the 


kinematic  variables  iT;  ,  and  w'*'  may  be  derived.  Elimination  of  from  0  using 
(172)  results  in 
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(173)  to  eliminate  <(>^1^  from 
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where  and  <f>^'  need  not  be  differentiable.  If  we  eliminate  w'^’  from  (183),  the 
extended  formulation  w’hich  does  not  require  continuous  differentiability  in  any  of  the 
kinematic  variables  is  realized. 
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In  the  context  of  use  in  the  finite  element  procedure,  fl,  leads  to  displacement 
formulations  while  leads  to  stress  formulations.  Evidently,  other  extensions  of  the 
general  variational  principle  than  the  ones  presented  here  are  possible.  For  example, 
elimination  of  derivatives  of  certain  force  resultants  from  fl,  and  fl,  or  derivatives  of 
kinematic  variables  from  fl^  and  results  in  so-called  mixed  formulations. 


61 


4.3.4  Some  Specializations 

If  the  admissible  state  is  constrained  to  satisfy  some  field  equations  and/or 
boundary  conditions,  certain  specialized  forms  of  the  variational  principle  are  realized. 
This  procedure  is  used  to  reduce  the  number  of  free  variables  in  the  governing 
function.  Also,  certain  assumptions  in  the  spatial  or  temporal  variation  of  some 
variables  can  lead  to  approximate  theories.  Some  specializations  of  the  extended 
variational  principles  stated  in  the  previous  section  are  presented  below. 

Specialization  of  the  function  flj,  to  the  case  where  (146)  and  (169)-(171)  are 
identically  satisfied  i.e.  if  displacement  w'  is  constant  through  the  thickness  of  plate 
and  the  iump  discontinuities  in  'displacement'  components  over  internal  surfaces  are 
identically  satisfied  leads  to 
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Since  w  uas  assumed  to  be  constant  over  the  tliickness  of  layer,  its  continuity  in  the 
interface  implies  that  the  lateral  displacement  is  constant  through  the  thickness  of 
plate.  If  we  further  specialize  O,  to  identically  satisfy  the  kinematic  relations 
(125)-(128), 
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Furtherfore,  if  the  interface  continuity  condition  of  the  in-plane  displacements  (145) 
and  the  displacement  boundary  conditions  (163)-(165)  are  identically  satisfied,  fig  in 
(186)  reduces  to  the  potential  energy  type  variational  principle.  Here,  two  different 
forms  of  expression  are  possible,  depending  on  w’hich  of  the  variables  is  eliminated. 
When  is  eliminated,  it  becomes 
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\'.f)ith  IS  the  vuriational  principle  for  Sun's  [l973]  theory.  On  the  other  hand, 
eliminating  </>'*',  we  have 
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(188) 


which  is  the  variational  principle  for  Srinivas  (197.^)  and  Seide's  [1980]  theories.  In 
connection  with  finite  element  formulation,  it  is  worth  noting  that  use  of  n,o  is  more 
convenient  if  in-plane  streching  of  individual  layer  needs  to  be  specified. 

Clearly,  a  large  number  of  other  specializations  from  other  extended  variational 
principles  are  possible.  We  desist  from  an  attempt  to  make  a  comprehensive  catalogue 
of  such  possibilities. 


4.4  COMriJiMhNTARY  VARIATIONAL  PRINCIPLES 

4.4.1  General 

An  alternative  procedure  to  set  up  variational  principles  governing  the  problem  is 
to  write  the  operator  equations  in  complementary  form  instead  of  the  direct 
formulation  (147).  In  this  formulation,  it  is  assumed  the  kinematic  relations  are 
satisfied.  Al-Ghothani  [1986]  presented  the  complementary  formulation  of  laminated 
composite  plate  for  the  dynamic  case  using  a  discrete  laminate  theory  and  discussed 
various  specializations  of  the  extended  variational  principles.  In  this  section,  we  present 
the  complementary  form  of  the  field  equations  given  in  the  previous  section.  Except 
the  constitutive  equations  for  transverse  shear,  the  formulation  is  the  same  as  the  one 
given  by  Al-Ghothani  [1986].  Since  an  extensive  di'^cussion  on  the  extended  principles 
and  various  specializations,  some  of  which  led  to  the  variational  principles  of  various 
approximate  theories,  has  been  given  in  [Al-Ghothani  1986],  those  investigations  will 
not  be  repeated  here.  However,  some  extensions  of  the  general  complementary 
variational  principle  and  specializations  which  are  not  included  in  [Al-Ghothani  1986], 
but  are  interesting  are  presented. 
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4.4.2  Complementary  Form  of  the  Field  Equations 

Assuming  the  kinematic  relations  (125M128)  are  satisfied,  the  field  equations 
(129)-(131),  (137)-<139)  and  (145),(146)  may  be  written  in  the  self-adjoint  matrix  form 


t*  t*t^  0  0  0  0 
0  0  /*  0  0  0 

-t*  0  0  0  0  0 
0  0  -t*  0  0  0 


(189) 


(190) 

(191) 

(192) 
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E  = 

‘<j 


0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 


0  0  0  0  0 


t.j=  1.  2. 


n 


(193) 


u,  =  [w^  ,  ,  w  ,  .  Q,  J.  *=  1.  2 . n 

=!  =  [  ■/■“’.  7'“']  fc=  1.2, . 

/  =  +  /♦G  +  0,0,0] 

4  i\  n  I"*  3 

==  [<*7’1"'.  0,  /*7’3"',  0,0,0] 


^  [j,  Q  Q  Qj 

f?  o’  no’  3 


n,  =  8  c-^  +  8  A 
'  97  “^98 

n  =  8  -4-  +  S« 

^  “>  d0  da 


n-l 


(194) 

(195) 

(196) 

(197) 

(198) 

(199) 

(200) 


The  operator  matrix  in  (189)  is  self-adjoint  in  the  sense  of  (A-25)  if  the  bilinear 
mapping  defined  in  (159)  is  used.  This  self-adjoint  form  of  the  field  equations  (189)  is 
different  from  Al-Ghothani's  [1986]  in  that  it  includes  matrices  E,^  representing  the 
coupling  of  transverse  shear  constitutive  relations  between  layers  based  on  the  consistent 
shear  theory  developed  in  Section  III. 


4.4.3  General  Complementary  Variational  Principle 

For  the  operator  matrix  equation  (189),  the  governing  function  is  defined, 
following  (A-29),  as 


n-  I 


j  = 


<U  ,  A,U,  > 


<  u 


r  „  - 


+  £  <uj,c=,  ,> 
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n  n 


i-l  y-1  it^l 

n-l  n- 1 

+  £  <=:,  b[u.>  ,„  +  £  <sl.<fu,..  > 


,(*) 


i-1 


*-l 


n 

+  2Z  << .  r,  >  j.,  -  2<  <•  S„  >^„  +  2<  U".  5,  > 


,(;i) 


/  -  1 


■^Boundary  Terms  +  Internal  Jump  Terms  (201) 

Substituting  (190)-(2(K))  into  (201),  the  explicit  form  of  the  governing  function  is 
obtained. 
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(202) 


As  in  direct  formulation,  it  can  be  shown  [Al-Ghothani  1986]  that  the  Gateaux 
differential  of  J  vanishes  if  and  only  if  the  field  equations  (189)  along  with  the 
boundary  conditions  and  the  internal  jump  discontinuity  conditions  are  satisfied. 


4.4.4  Extended  Complementary  Variational  Principles 

Following  the  principles  and  methodology  presented  previously,  it  is  possible  to 
de\elop  extended  variational  principles  for  the  complementary  form  of  the  field 

equations  (189)  as  well.  Relations  (172)-(174)  can  be  used  to  eliminate  some  of  the 
operators  from  (202).  Simultaneous  use  of  those  relations  to  eliminate 

and  results  in 
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Here,  the  force  resultants  need  not  be  differentiable.  Alternatively,  eliminating 
^*3’ *^0.3  from  (202), 
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(204) 


-2<w'\ -ig!,!))  >k  -  2<V-^,  >  ,,, 

^S/  ^2i 

-2<m'^^,  -  2<Qf,  >^J 


In  addition  to  extensions  illustrated  above,  it  is  obvious  that  numerous  other  extended 
principles  are  possible  [Al-Ghothani  1986]. 


4.4.5  Some  Specializations 

As  in  the  direct  formulation,  some  specializations  of  the  complementary  extended 
variational  principles  are  possible  by  requiring  that  certain  field  equations  and  or 
boundary  conditions  be  identically  satisfied. 

If  we  assume  that  the  equilibrium  equations  (129)-(131)  are  identically  satisfied, 
results  in 
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If  we  further  specialize  (205)  to  satisfy  the  stress  boundary  conditions  on 
.S’i‘'.-^3*’  and  5*/’,  and  the  internal  jump  discontinuity  conditions 

J.  =  £<  <“"•  u,  +  «f>V^  u,  +  ,,  +  2<u[[\ , 


,0)  (X), 


+  £  I  <cr. 


»(*) 


i-1  yl 
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4.5  DlSCZfSSION 

Based  on  the  discrete  laminated  plate  theory,  which  accounts  for  the  effect  of 
transverse  shear  deformation  in  a  consistent  manner,  a  systematic  development  of 

variational  principles  for  dynamics  of  linear  elastic  composite  laminated  plate  has  been 
presented.  The  direct  as  well  as  complementary  formulation  are  considered. 

Complementary  self-adjoint  form  of  the  field  equations  is  the  same  as  the  one 
presented  by  Al-Ghothani  [1986],  except  for  the  coupling  terms  of  transverse  shear 
constitutive  equations  between  layers  which  have  been  introduced  in  the  consistent 
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shear  JeCormable  theory  presented  in  Section  111.  Nonhomogeneous  boundary  conditions 
and  internal  jump  discontinuities  have  been  explicitly  included  in  general  variational 
principles.  Allowance  of  jump  discontinuity  terms  in  variational  formulation  is 
necessary  in  the  context  of  direct  approximation  in  finite  element  spaces  since  the  space 
of  approximants  may  not  be  sufficiently  smooth.  Also,  extensions  of  the  general 
variational  principles  through  elimination  of  certain  field  operators  and  specializations 
by  restricting  some  of  the  field  equations  and/or  boundary  conditions  to  be  identically 
satislied  have  been  proposed,  figure  2  and  figure  ^  diagrammatically  depict  possible 
extensions  of  the  general  variational  principles  based  on  the  direct  and  the 
complementary  formulation,  respectively.  In  either  case,  the  specializations  listed  in  this 
section  are  shown.  These  formulations  can  provide  a  basis  for  development  of 
alternative  approaches  to  approximate  solution  of  the  problem  and  also  for  development 
of  approximation  theories.  Evidently,  other  extended  forms  could  be  used  as  starting 
points  for  specialization. 
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:  Governing  Function 
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li^ure  2:  Family  of  Variational  Principles  by  Direct  Formulation 


9 


Figure  3:  F’amily  of  Complementary  Variational  Principles 


Section  V 


FINITE  ELEMENT  FORMULATION  OF  A  SPEOAL  DISCRETE 
LAMINATE  THEORY  OF  PLATES 


5.1  I  NTROIXICnON 

In  tlie  I'mite  element  priKcdure,  the  region  under  consideration  is  subdi\ ided  into  a 
t  inite  number  ol'  dis^unt  subregions  (elements),  and  the  field  variables  of  the  problem 
are  approximated  b\  functions  vchith  are  continuous  along  the  boundary  of  elements, 
but  have  limited  smcxtthness.  Consider  the  open  connected  region  R  in  an  Euclidean 
space  discreti?ed  by  a  finite  number  of  elements  (R,,  c=l,  .  .  m)  such  that 

R  =  Urn  M  R,  (207) 

1 

in  which  the  elements  satisfy  the  property 

R  =  if  e^f  (208) 

and  are  connected  at  a  finite  num(>er  of  nodal  points.  Here,  R  and  A,  denote  the 
closures  (if  R  and  R,  respectively. 

Since  the  field  variables  of  the  problem  are  represented  by  functions  which  may 
not  be  sufficiently  smtxith  over  the  region  R,  the  variational  principles  derived  in  the 
previous  section  may  not  be  valid  over  the  region  R.  However,  if  approximate 
functions  have  adequate  smoothness  over  each  element,  and  internal  discontinuities 
across  element  boundaries  are  explicitly  included,  they  are  valid.  For  such  case,  we 
define  the  governing  function  over  the  region  R  as 


80 


(209) 


« =  z  «< 

<  1 

where  Q,  (e=l,2,....jn)  is  the  set  of  functions  governing  the  problem  over  indicated 
element  R^.  Sandhu  [1976]  showed  that  Gateaux  differential  of  Q  in  (209)  vanishes  if 
and  only  if  the  field  equations,  the  boundary  conditions  and  the  continuity  conditions 
across  the  interelement  boundaries  are  satisfied.  If  there  are  discontinities  across  the 
interelement  boundaries,  the  actual  jump  quantities  need  to  be  explicitly  included  in  fl, 
[S<.ndhu  1976], 

In  the  following,  we  present  a  finite  element  formulation  based  on  the  variational 
principle  using  the  governing  function  ft,  which  is  defined  in  terms  of  displacement 

field  variables  and  w,  and  is  a  specialization  of  the  extension  of  the  general 

variational  principle  to  cases  where  M^g,  may  not  be  continuously 

differentiable.  The  specialization  is  to  the  case  where  >v'*’  =  w  for  all  k,  the 
strain-displacement  equations  are  identically  satisfied,  the  in-plane  displacements  are 
ccmtinuous  and  the  specified  displacement  boundary  conditions  are  identically  satisfied. 


.^.2  SPATIAL  DISCRiaiZATION  OF  GOVERNING  FVNCl'ION 

The  governing  function  ft,,  satisfying  the  kinematic  relations  and  displacement 
continuities  in  the  laminar  interfaces,  can  be  rewritten  as  the  sum  of  functions  ft^ 

(e=1,2 . m)  defined  on  each  element  of  the  finite  element  representation  as 
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where  ^u\  and  arc  the  intersections  of  the  boundary  surface  of  an  element 
with  S^,*\  S^j*'  and  respectively.  For  spatial  discretization  of  this  function,  it  is 
assumed  that  the  field  variables  are  interpolated  over  an  element  domain  as 

u'‘’(x.t)  =  H^(x)U(0  (211) 
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(212) 


w(x,i)  =  H^(x)W(t)  (213) 

W 

where 


and  U(«),  <!>''’(/ ),  W(t)  are  the  vector  functions  of  time  defined  at  the  nodal  points  and 
are  respectively,  the  matrices  of  spatial  interpolation  functions  for  the  field 
variables  indicated  by  s  bscripts.  Also,  the  generalized  strains  may  be  expressed  as 
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and  T^.  T^,  and  T,  are,  respectively,  the  transformation  matrices  derived  from  the 
displacement  interpolation  functions  and  H„.  by  suitable  differentiation  and 

reorganization  of  terms. 

Substituting  (21l)-(216)  into  (210),  the  spatially  discretized  governing  function  is 
obtained. 
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<-  I  XI 


=  -  U  M  *  U  +  2U  +  Yt  *yt<t>' 
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09/  WW  U0^ 


+2Tr(<I)'/M'^-*<l>^*'  +  +  2f*uV*j/  Tf  <I>' 

i  00^  ukA  U0.A  / 


/  I  /-I 


'  00'\  J  00/'  //0/ 


+  2/*  V/C^D'/kY'  *<t>"'} 

4L^  I  00/> 


/  .  r/K<‘ V  aO)\ 


X • 1  y- 1 


X'l  i»l  !•! 

+/‘(<I>'*y*g^’  +  (4>^^Y*b^*'  +  +  W^*b‘*^ 

y  3  .. 

n 

+2rM-U^*T'’  +  W^*(p''^-p‘')  +  U^*t"’  +  } 


-It*  Vi  u'  ♦r"’  +  y  J  (<I)'/‘r'‘,’  +  («I>“Y*r‘'’  +  W^*r'^^}  (217) 

'  I  /  n2  m  y 


here 


»l‘'=  f  H  dR 

M".  =  r  H  P'"  'h^  xfR 

Utpp  J  L  0 

?- 

f 

m"’  -  farvjf 

I'W  I  't>  V 

:fx> 


(218) 


(219) 


(220) 
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m'“  = 


— ■ 


f  H 

I  u  <t> 

;?*> 
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f 

f 

/H  dR 

“ 


dR 


dR 


m"  = 


k'>  r 

uuA  J 


T  V'V  dR 


k'*’  =  r  t 
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T  dR 
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T  dR 
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'‘<t^A 


I 


T  dR 
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f, 
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f 
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k  I  5  S 

r  T  </j? 

?*) 


.SW0 


(221) 

(222) 

(223) 

(224) 

(225) 

(226) 

(227) 

(228) 

(229) 

(230) 

(231) 

(232) 
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k‘'^=  r  H 

'W  10  <P 


f'‘'  =  r  H  d/? 
■ 

fj  ’  =  r  H^Fj  ’  rfF 
3*) 


g"  =  H^G 


h‘  =  /  H  F''^//e 


^'. '  =  / 

J*) 

b<*>=  f 

it) 


x'"’ 


.«) 


b  '  =  ^  H  X  dR 


At) 


b'  '  =  I  H^Y  '  </F 
0 


.U) 


/“ 

=  /«. 

?») 


«) 


,(i) 


b'  '  =  I  H  Z 


(f’,T''’)=  r  H.  (T".  t'"') 


•f " 


Jn) 


r.('l) 
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(p”.p"^=  f  lljTl,fj)dR 

jr*’ 

=  r  H  n*  \ii 
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(233) 

(234) 

(235) 

(236) 

(237) 

(238) 

(239) 

(240) 

(241) 

(242) 

(243) 

(244) 

(245) 
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= 

u2 


f 

A)  ^ 


ds 


(k)  j 

m  ds 


=  f  ^ 

3*) 


rU  )  ^ 
9 


and 


(246) 


(247) 


(248) 


X«> 

=  [x“’, 

,  Y 

yUY 

p(i) 

= 

g'*’ 

g[^Y 

T“  = 

=  [r^ .  r°/ 

=  It\"\ 

yny 

= 

=  [< 

\  Xf';’ 

(< ) 

q 

=  [0‘;\ 

0*;’] 

(249) 

(250) 

(251) 

(252) 

(253) 

(254) 

(255) 

(256) 

(257) 


In  the  above  expressions,  A'*’,  B“’  and  denote,  respectively,  the  matrices  of  the 
stiffness  aI*J^ 

*  ^a0y6  and  For  convenience,  we  rewrite  the  quantities  involving 

summations  in  (217)  as 

£{2U"Ml‘^/^r<I.'}  =  2U^£(  (258) 

<  =  1  1=1  <-i  i-k*t 
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n  X  - 1 
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n~  1 
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n  X-1 


fi-l 


£(2£.,(*yK"  •♦<*>)  =  2£i.(*"VM  £k;„*') 

t-l  i-l  *-l  i-i+l 

Using  above  equalities,  the  discretized  function  (217)  becomes 

=  -<  ZuX*>U  +  ZKi.*,>****‘*^  + 


(263) 
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T^(k)  *  j^(^) 


n-  1 


kW>J, 


t-l 


i-l 


i-l 


i-i-t  I 


-t*  Z^’W^Z*^lvv*  ^ 

*-i  j-i  j-i  j-i 


*-i 
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+2t*{-U^*T°  +  W^»(p''’-p‘’)  +  +  ^//«I)'/*t'"'  } 

1-1 


(264) 


*-i 


1-1 


5J  SEMIDISCRETE  EQUATIONS  OF  MOTION 

To  obtain  the  semidiscrete  equations  of  motion  of  laminate  plate,  it  is  convenient 
to  rewrite  the  spatially  discretized  variational  principle,  (264)  in  matrix  form  as 

n  =-X^S*X  +2X^»R  (265) 

€  e  e  e  e  e 

where 
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S,.  1  5,,  1  0 
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K, '  I 
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0  I  5^31533 


(266) 


Here,  elements  of  the  submatrices  of  S,  and  the  load  vector  are,  explicitly. 


S,  ,  =  + 

1 1  uu  utiA 


Z  +  .  j  =1.  2.  •  •> 

*-y+i 
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l-^l 

The  spatially  discretized  governing  function  for  the  global  system  is  given  by 

m 

n  =  = -X^S»X  +  2X^*R  (26?) 

o- 1 

where  X  is  the  vector  of  values  of  field  variables  at  the  system  nodal  points,  R  is  the 
set  of  corresponding  'forcing'  quantities  and  S  is  the  system  matrix  corresponding  to  S, 
for  an  element.  Here,  summation  in  (267)  indicates  matrix  assembly  following  usual 
(direct  stiffness)  procedure.  Vanishing  of  the  differential  of  n  in  (267)  with  respect 
to  \  gives  the  set  of  equations: 

S\=M\  +  /*KX  =  R  (268) 


where 
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The  semidiscretized  equations  of  motion  for  the  entire  system  can  be  obtained  by 
differencing  (268)  twice  with  respect  to  time. 

MX  +  KX  =  R  (269) 

where  mass  matrix  M  and  stiffness  matrix  K  are  symmetric.  Furthermore,  M  and  K 
are  positive  definite  and  semi-positive  definite,  respectively. 

5.4  I'REE  VIBRATION  ANALYSIS 

For  free  vibration  analysis,  the  load  vector  R  is  set  to  zero  and  (269)  becomes 

MX  +  KX  =  0  (270) 

Assuming  harmonic  motion  of  the  system,  i.e.,  assuming  the  solution  by  X  = 
where  ^  is  the  amplitude  vector,  w  is  the  natural  frequency  and  i=sfA,  we  obtain 
generalized  eigenvalue  problem 

(K-«‘^^M)«^-  =  0  (271) 

Before  applying  the  boundary  constraints,  this  eigenvalue  problem  has  three  zero 
eigenvalues  corresponding  to  the  rigid  body  motions.  If  equations  corresponding  to  the 
constrained  nodes  are  removed  before  solving  the  eigenvalue  problem,  matrix  K  becomes 
positive-definite  and  consequently  all  eigenvalues  are  positive  real  and  corresponding 
eigenvectors  representing  vibration  modes  are  orthogonal  with  respect  to  the  mass 
matrix  M. 
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5J  SPECIALIZAI  ION  TO  THE  STAl'IC  CASE 

A  specialization  to  the  static  problem  can  be  made  by  dropping  the  inertial  terms 
from  the  set  of  equations  (269).  The  resulting  system  of  linear  algebraic  equations  is 
KX  =  R  (272) 

Equation  (272)  can  be  modified  to  account  for  known  boundary  constraints  and 
subsequently  solved  for  the  unknown  vector  X  comprising  kinematic  quantities. 

5.6  EREE  EDGE  DELAM I N Al'ION  SPECIMENS 

For  the  case  of  a  free-edge  delamination  specimen,  the  stresses  at  the  free-edge 
are  known  to  be  zero.  As  the  field  variables  in  the  present  formulation  are  the  nodal 
values  of  displacements  and  rotations,  these  stress-free  constraints  cannot  be  directly 
applied  to  the  system  (272).  Since  in  this  formulation  stresses  are  secondary  variables 
to  be  computed  from  the  obtained  displacements,  appropriate  refinement  of  the  mesh 
near  the  edge  boundary  would  be  needed  to  realize  as  close  to  zero  as  possible. 
However,  as  discussed  in  Section  3.7.4,  different  constitutive  equations  (122)  apply  at 
the  free  edge  for  the  case  of  transverse  shear.  This  leads  to  two  sets  of  constitutive 
relations,  one  applicable  in  the  interior  of  the  laminate  and  the  other  at  the  free  edge. 
Hence,  in  evaluating  the  integrals  in  (231)  through  (233),  (108)  were  used  for  elements 
lying  entirely  within  the  interior  of  the  laminate.  For  elements  having  one  edge 
coinciding  with  the  free  edge,  values  of  at  each  of  its  nodes  were  chosen  from 

either  (108)  or  (122)  depending  on  whether  the  node  was  on  the  free  edge  or  in  the 
interior.  Then  the  value  of  A'*'*  at  the  gauss  points  could  be  determined  by 
interpolating  the  nodal  values  using  some  interpolating  functions.  For  the  examples  of 
application  given  in  this  report,  the  same  interpolation  was  used  as  adopted  for  the 
displacements. 
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5.7  CAIjCULATION  OF  S'FRFJSSES 


After  obtaining  the  displacement  solution,  secondary  computations  to  obtain  the 
stress  components  can  be  carried  out  in  the  following  manner. 

S.7.1  In-PIane  Stress  Components 

Using  the  known  nodal  displacements  in  conjunction  with  the  interpolating 
functions,  strains  can  be  obtained  at  any  desired  location  over  an  element.  Subsequently 
using  (25),  and  noting  that  e%'  =  0  in  this  formulation,  the  components  of  in-plane 
stress  can  be  calculated.  Alternatively,  (28)  may  be  used  to  obtain  the  stress  resultants 
and  (75)  used  to  obtain 


5.7.2  Transverse  Shear  Stresses 

Having  the  solution  for  the  nodal  displacements  and  rotations,  the  resulunts  Q'J'’ 
can  be  determined  from  the  constitutive  relations  (108).  Subsequently,  use  of  (92)  will 
yeild  Xh  hence  The  transverse  shear  stresses  can  then  be  evaluated  from 
(81). 

For  the  case  of  the  free-edge  delamination  specimen,  interpolation  over  the  edge 
elements  as  described  in  Section  5.6  was  done  to  obtain  and  7^^*'  at  the  desired 
locations  over  the  element. 


5.7.3  Transverse  Normal  Stress 

To  compute  the  normal  stress  cr„,  the  equilibrium  equation  (77)  is  integrated  w.r.t 
Xj  to  obtain 


(273) 


using  (81 ), 
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(274) 


_  r  [  ^  ‘  V‘*  ^  ]  dx^‘ 

33  3  I  ^1  *'o,a  ^2  o.o  ^3  oi.a  3 

^  0 

Substituting  for  and  {3  from  (81)  and  evaluating  the  integrals  leads  to 
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t. 
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a.o 


(275) 
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Section  VI 


APPLICATIONS 


6.1  INTRODUCTION 

■|'he  formulation  presented  in  the  preceding  section  was  incorporated  in  a  computer 
program.  A  nine-nixied  Heterosis  element,  shown  in  I'igure  H.l,  was  used  for  finite 
element  analysis.  For  verification  of  the  code,  the  program  was  used  to  solve  the 
problem  of  a  three-layer,  simply  supported,  square  sandwich  plate  with 
isotropic/orthotropic  outer  skin  layers.  The  results  were  compared  with  an  exact  series 
solution.  TTiis  example  is  the  same  as  used  by  Mawenya  [1974].  The  formulation  was 
then  applied  to  the  solution  of  a  Free  Edge  Delamination  (FED)  problem.  A  four-layer 
coupon  under  uniform  axial  strain  whose  solution  has  been  presented  by  Pagano{l978] 
was  considered.  Displacements  and  stresses  along  the  midsection  were  computed  and 
compared  with  Pagano's  results. 

Though  the  constitutive  relations  for  transverse  shear  have  been  derived  using  the 
equations  of  static  equilibrium,  the  problem  of  free  vibration  of  a  simply  supported 
laminated  plate  was  also  solved  by  including  inertial  terms  in  the  finite  element 
formulation.  Developing  constitutive  equations  for  the  coupled  theory,  allowing  for  the 
inertia  terms,  Schoeppner  [1990]  has  shown  that  the  effect  of  inertia  on  constitutive 
relations  decreases  with  decrease  in  layer  thickness.  Thus,  for  sufficiently  small  layer 
thickness  (any  lamina  can  be  arbitrarily  replaced  by  a  suitable  number  of  sublayers) 
the  constitutive  relations  developed  in  Section  IV  would  be  applicable.  The  natural 
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frequencies  of  free  vibration  of  a  sandwich  plate  were  computed  and  the  fundamental 
frequency  compared  with  an  exact  elasticity  solution  by  Srinivas  [1970]. 


6.2  PLATE  ANALYSIS 

A  three-layer,  simply  supported,  square  plate  uniformly  loaded  in  the  transverse 
direction  was  considered  for  analysis.  Two  separate  cases  in  which  the  outer  layers 
were  respectively  isotropic  and  orthotropic  were  considered.  The  plate  dimensions  and 
material  properties  were  the  same  as  in  a  similar  example  considered  by  Ma\xenya 
[1974],  and  were 

Plate  Dimensions 

Length  of  each  side  =10  in. 

Thickness  of  outer  layers  =  =  0.028  in. 

Thickness  of  core  t^^^  =  0.75  in. 

Material  Properties 

a).  Isotropic  Outer  Layers 

£•"’  =  =  10''  ib/in.^ 

=  0.3 

g' ^’ =  3x  lO'*  a»/in.^ 

•  b).  Orthotropic  Outer  Layers 
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I-\'l  =  /i',’’  =  lo'  Ib/in.^ 

^22  =  =  ■<  ><  10*“ 

g\'1  =  =  1.875  X  lO**  Ib/in.^ 

v\'l  =  =  0.3 

g',’,'  =  3  X  lo"  Ib/itu 
g\'^  =  1.2x10'*  Lb/in.' 

Lateral  Loading 

^  =  1.0  Ib/in.^ 

Due  to  the  symmetry  of  the  problem,  only  one  quadrant  of  the  plate  needed  to  be 

considered  for  finite  element  analyses  under  static  conditions.  A  typical  discretized 
quarter  along  with  associated  boundary  conditions  is  shown  in  Figure  4.  Discretization 
was  done  using  1x1,  2x2,  4x4  and  8x8  meshes. 

The  results  for  maximum  lateral  deflection  at  the  plate  midpoint  are  shown  in 
Table  5  and  Table  6  along  with  the  exact  series  solutions  [Mawenya  1974]  and  a 

comparison  with  results  obtained  by  Moazzami  et  al.  [l99l]  using  the  conventional 
Discrete  Layer  theory.  The  CPU  time  used  on  a  Cray  Y-V1P8/864  is  also  listed.  It 

shou'd  be  noted  here  that  the  Heterosis  element  used  is  nine-noded  except  for  the 
lateral  displacement  degrees  of  freedom,  for  which  it  is  eight-noded.  Moazzami's  results, 
on  the  other  hand,  were  obtained  using  a  fully  nine-noded  element. 

These  results  were  obtained  using  a  generalized  plane  stress  assumption  for  the 
purpose  of  matching  them  with  the  series  solution.  Solving  (27)  for  e,,’  and 

substituting  into  (25)  and  (27), 
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v=4)f =0 


0/2 


u  =  =  w  =  0 


« 


Hgure  4.  Discretized  Quadrant  of  a  Simply  Supported  Plate 
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Table  5 


Maximum  Lateral  Deflections  for  the  Isotropic  Case 


No.  Of 
Elements 

CPU  (Secs. ) 

Maxm .  Def 1 . 
from  Code 

Max.  Defl. 
Moazzaml  et  al . 

Ixl 

0.113 

0.00069822 

0.00073539 

2x2 

0.411 

0.00073456 

0.00074087 

4x4 

1.733 

0.00073537 

0.00074008 

8x8 

9.912 

0.00073538 

0.00074006 

Series  Soln. 

[Mawenya  1974] 

0.00074000 
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Table  6 


Maximum  Lateral  Deflections  for  the  Orthotropic  Case 


No.  Of 
Elements 

CPU  (Secs. ) 

Maxm.  Defl. 
from  Code 

Max.  Defl. 
Moazzami  et  al. 

1x1 

0.114 

0.0011452 

0.0012140 

2x2 

0.414 

0.0012063 

0.0012224 

4x4 

1 . 783 

0.0012074 

0.0012216 

8x8 

9.909 

0.0012075 

0.0012216 

Series  Soln. 

[Mawenya  1974] 

0.0012300 
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Assuming  generalized  plane  stress  state  in  a  lamina,  i.e., 


0 


(279) 


the  constitutive  equations  of  bending  and  stretching  are  obtained  as  in  (28),  where  now 


(280) 


We  note  here  that  the  assumption  of  Uj*’  consunt  over  the  thickness  of  each  layer 
implies  ej3’  =  0.  We  also  note  that  (25)  and  (27)  would  not  contain  and  hence  there 
would  be  no  need  for  its  elimination.  The  in-plane  constitutive  equations  can  be 
obtained  as  was  done  in  (28)  without  any  assumptions  on  However,  for 

comparison  with  the  series  solution  based  on  Kirchoffs  theory,  (280)  were  used.  The 
choice  of  the  two  constitutive  relations  is  easily  incorporated  into  the  program  with 
the  inclusion  of  a  few  lines  of  code. 

The  assumption  of  generalized  plane  stress  state  defined  by  (279)  was  first  used 
by  Mindlin  [1951]  for  homogeneous  plate  theory  and  later  by  Whitney  [1970]  for 
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laminated  plate  theory.  They,  however,  assumed  the  integral  of  over  the  entire 

thickness  of  the  laminate  to  vanish  as  opposed  to  over  each  layer  indicated  by  (279). 
Hong's  [1988]  assumption  (279)  is  clearly  invalid  for  a  multi-layered  laminate. 
However,  for  a  single  layer  plate  or  the  sandwich  plate  with  thin  outer  layers  under 
consideration,  it  is  acceptable. 


6J  APPLICATION  TO  FREE  EDGE  DELAMl N AJ'JON  {FED)  SPECIMEN 

Pagano[l978]  presented  solutions  for  the  free-edge  problem  of  a  coupon  subjected  to 
uniform  axial  strain  e.  A  four-layer  laminate  as  shown  in  Figure  5  was  considered. 
The  material  properties  were: 

£, ,  =  20.0  X  lO**  ,  £22  =  ^33  =  ^ 

Gj2  =  G,3  =  G23  =  0.85X10^  {psi)  (281) 


*'.2  =  ‘'13  =  *^23  =  0-21 


To  match  the  ratios  of  the  spatial  dimensions  with  those  used  by  Pagano,  the  laminate 
length  :  width  :  layer  thickness  (i.e.,  L  :  2b  :  h)  was  selected  as  11  :  3  :  0.1875.  Two 
stacking  sequences  viz.,  [-^45/-45],  and  [0/90],  were  considered. 


6.3.1  Angle  Ply  Specimen 

In  the  X  direction  the  specimen  was  discretized  using  11  equal  sized  elements.  In 
the  y  direction  a  sequence  of  mesh  refinement  was  done  by  starting  with  three  equal 
elements  and  succesively  subdividing  the  elements  along  the  two  traction-free  edges 
further  into  three  equal  elements.  Using  this  pattern,  3x11,  7x11,  11x11  and  15x11 
element  meshes  were  used.  This  scheme  is  illustrated  in  Figure  6  where  refinement 
from  the  3x11  to  the  7x11  mesh  is  depicted.  A  nine-noded  Heterosis  element  shown 
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in  Figure  B.l  was  used  for  the  finite  element  analysis. 

The  results  for  the  inplane  stresses  cr„  and  for  the  various  meshes  are 

compared  with  Pagano's  solution  in  Figure  7  and  Figure  8.  These  stresses  are  plotted 
across  the  middle  surface  of  the  top  layer  for  comparison  with  Pagano's  solution.  The 
stresses  are  shown  over  half  the  width  of  the  coupon,  where  (y-fc)/t  =  0,  1, 
respectively  represent  the  center  and  the  free  edge  of  the  laminate. 

The  result  for  the  3x11  mesh  shows  a  jump  across  the  element  boundary.  It 
should  be  noted  that  as  only  half  the  laminate  width  is  represented  in  the  plots,  the 


3x11  mesh  is  represented  by  1-j  elements.  For  the  chosen  element,  in-plane  strains,  and 


consequently  the  stresses,  are  linear  over  the  element.  The  edge  element  in  the  3x11 
mesh  attempts  to  model  both  the  sharply  varying  stresses  near  the  free  edge  and  the 
relatively  flatter  ones  near  the  center  by  a  single  linear  fit.  Hence  the  stresses  near  the 
inner  boundary  of  the  edge  element  tend  to  be  overestimated.  This  explains  the  kink 
in  the  3x11  solution  across  the  element  boundary.  With  mesh  refinement  the  kink  was 
observed  to  reduce  for  the  7x11  mesh  and  was  not  noticeable  in  the  solution  for  the 
11x11  clement  mesh.  Moreover,  it  was  observed  that  further  refinement  to  a  15x11 
mesh  did  not  show  "substantial”  improvement  in  the  solution.  A  comparison  of  the 
total  CPU!  time  for  the  different  meshes  using  a  CRAY  X-MP/28  is  given  in  Table  7. 

The  results  in  Figure  7  and  Figure  8  show  that  the  finite  element  results 
overestimate  Pagano's  solution  for  O',,  and  by  about  4  percent  at  the  center  and  12 

percent  at  the  free  edge.  The  same  amount  of  error  is  observed  at  the  free-edge  in 
the  longitudinal  displacements  plotted  in  Figure  9  for  the  "optimum”  11x11  mesh. 
The  transverse  shear  stress  a,,  computed  using  (81)  is  shown  in  Figure  10. 
Comparison  with  Pagano's  exact  solution  shows  that  the  numerical  results  grossly 
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Figure  7;  Distribution  of  X-stress  along  Center  of  Top  Layer  (Angle-Ply). 
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Table  7 


Comparison  of  CPU  time  for  different  Finite  Element  Meshes. 


F.E.M  Mesh 

Total  CPU  time  on  a  CRAY  X-MP/28 

(Seconds . ) 

3x11 

7.798 

7x11 

23.424 

11x11 

42.664 

15x11 

68.921 

underestimate  the  shearing  stress. 

To  further  enhance  the  accuracy  of  the  solution,  refinement  in  the  thickness  direc¬ 
tion,  i.e.  the  division  of  each  layer  into  sublayers,  was  attempted.  Each  layer  was 
divided  into  three  sublayers  of  thickness  h/3.  Further  subdivision  into  five  layers  was 
done  by  dividing  the  sublayer  of  interest  (i.e.,  the  one  containing  the  location  at  which 
stresses  were  desired)  into  three  equal  layers  of  thickness  h/9.  To  save  on  comj/Jta- 
tional  effort,  the  symmetry  of  the  problem  was  exploited  by  considering  only  the  top 
two  plies  and  specifying  the  transverse  displacements,  w,  at  the  middle  surface  of  the 
laminate  equal  to  zero.  Hence,  solutions  for  the  cases  of  N  =  2,  7V=6  and  N=10  were 
obtained,  where  N  represents  the  total  number  of  resulting  layers;  N~2  being  the  case 
with  no  sublayers.  The  scheme  for  sublayer  division  is  depicted  in  Figure  11. 

The  results  for  inplane  stresses  for  different  N  using  a  11x11  element  mesh  are 
compared  with  Pagano's  results  in  Figure  12  to  Figure  15.  A  comparison  of  CPU  time 
for  various  number  of  sublayers  is  given  in  Table  8.  Table  7  and  Table  8  show  that 
for  the  problem  using  the  11x11  mesh,  exploiting  the  symmetry  by  considering  only 
two  plies  reduced  the  CPU  time  from  42.664  to  15.852  seconds. 

Figure  12  and  Figure  13  show  significant  improvement  in  the  results  for  inplane 
stresses  at  the  free  edge  with  increasing  number  of  sublayers,  though  the  stresses  at 
the  center  of  the  laminate  remain  essentially  unchanged.  The  results  for  axial 
displacement  are  shown  in  Figure  14  and  show  that  the  displacements  obtained  using 
N=6  and  N=10  agree  closely  with  Pagano's  solution. 

The  transverse  shear  stress  tr,,  at  the  ±45"  interface  was  computed  using  (81)  and 
is  shown  in  Figure  15.  With  increasing  N,  the  computed  free-edge  stress  steadily 
increased  in  magnitude.  Pagano  [l978]  observed  that  a,,  at  the  edge  grows  with 
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Axial  Displaceme.ii  Across  Top  Surface  (Angle-Ply) 
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Table  8 


Comparison  of  CPU  Time  for  11x11  Mesh  with  Different  Num¬ 
ber  of  Sublayers. 


N 

Total  CPU  time  on  a  CRAY  X-MP/28 

(Seconds. ) 

N=2 

15.852 

N=6 

79.333 

N=10 

197.084 
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NODAL  DISPLACEMENTS 


Figure  14: 


Axial  Displacement  Across  Top  Surface  with  N=6,  10.  (Angle-Ply) 


increasing  N  and  could  not  determine  whether  a  finite  limit  was  approached  for  large 
N.  The  values  of  the  transverse  stress  close  to  the  free  edge  for  different  number  of 
sublayers  are  presented  in  Table  9. 

To  observe  the  effect  of  the  consistent  shear  treatment,  the  results  were  compared 
with  those  obtained  by  Moazzami  et  al.  [l99l]  using  the  discrete  laminate  theory 
along  with  a  constant  shear  correction  factor.  The  stresses  and  displacements  at 
element  centers  for  a  11x11  mesh  with  N=2  are  compared  in  Table  10  to  Table  13. 

Table  10  shows  that  the  results  for  <T,,  do  not  show  any  "substantial”  difference. 
The  results  for  cr,^,  a,,,  and  u^  in  Table  11  to  Table  13,  however,  indicate  that  ttiose 
obtained  from  the  present  theory  are  slightly  better  than  those  obtained  using  the 
existing  discrete  laminate  theory.  As  observed  in  Section  6.2,  only  marginal  improve¬ 
ment  was  expected  by  the  introduction  of  consistent  shear  coupling. 

In  summary,  it  was  observed  that  the  results  improved  to  a  certain  extent  only 
with  spatial  refinement  of  the  finite  element  mesh.  Convergence  after  the  11x11  mesh 
was  slow.  Further  enhancement  in  the  solution  was  obtained  by  dividing  each  ply  into 
sublayers.  As  the  inplane  displacements  are  assumed  linear  over  the  thickness  of  each 
layer,  division  into  sublayers  resulted  in  piecewise  linear  displacements  over  each  ply 
contributing  to  increased  accuracy. 
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Table  9 


Growth  of  XZ-slress  at 

Free  Edge  wUh  Increasing  N 

N 

The  Present  Study 

Pagano 

2 

0.416 

1.664 

6 

0.955 

2.213 

10 

1.450 

- 
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Table  10 


Comparison 

of  X -stresses  with 

those  Obtained  by  Moazzami  et  al. 

(y-b)/b 

Present 

Study 

Moazzami 

Pagano 

(N=6)* 

0.000 

3.0912 

3.0882 

2.902 

0.444 

3.0853 

3.0785 

2.90 

0.666 

3.0593 

3.0584 

2.89 

0.815 

2  9999 

3 . 0067 

2.78 

0.888 

2.9229 

2.9317 

2.67 

0.962 

2.7712 

2.7658 

2.25 

*  Note;  These  numbers  are  approximated  from  the  plots  in  Pagano  [197 
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Table  11 


Com  parison 

of  XY-stresses  with 

those  Obtained  by  Moazzarrd  et  al. 

(y-b)/b 

Present 

Study 

Moazzami 

Pagano 

(N=6)* 

0.000 

1 . 1793 

1.1995 

1  .  15 

0.444 

1 .  1668 

1 . 1921 

1  .  14 

0.666 

1 .  1137 

1 . 1526 

1.05 

0.815 

0.9923 

1 . 0504 

0.88 

0.888 

0.8348 

0.9023 

0.64 

0.962 

0.5248 

0.6409 

0.25 

*  Note:  These  numbers  are  approximated  from  the  plots  in  Pagano  [197 
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Table  12 


Comparison  of  XZ-stresses  with  those  Obtcdned  by  Moazzami  et  al. 


(y-b)/b 

Present 

Study 

Moazzami 

Pagano 

(N=6)* 

0.000 

-0.0030 

-0 

o 

o 

o 

1 

0.444 

-0.0143 

-0.0072 

-0.019 

0.666 

-0.0613 

-0.0387 

-0.09 

0.815 

-0.1752 

-0.1268 

CD 

CJ 

o 

1 

0.888 

-0.3046 

-0.2385 

CO 

to 

o 

1 

0.962 

-0.5465 

-0.4662 

-0.939 

*  Note:  These  numbers  are  approximated  from  the  plots  in  Pagano  tl97 
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Table  13 


Comparison  of  X -displacements  with  those  Obtained  by  Moazzami 

et  al. 


(y-b)/b 

Present 

Study 

Moazzami 

Pagano 

(N=6)* 

0.000 

-0 

-0 

0.00 

0.444 

-0.0143 

-0.0071 

-0.011 

0.666 

-0.0615 

-0.0387 

-0.055 

0.815 

-0. 1765 

-0. 1244 

-0.16 

0.888 

-0.3069 

-0.2338 

1 

o 

to 

05 

0.962 

-0.5507 

-0.4571 

-0.47 

*  Note:  These  numbers  are  approximated  from  the  plots  in  Pagano  [197 
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6.3.2  Cross-Ply  Specimen. 

The  second  example  of  the  FED  specimen  solved  was  a  [0/90]^  laminate.  The 

material  properties  and  laminate  geometery  were  the  same  as  before.  Results  for 
CTjj  and  <T„  were  computed  for  comparison  with  Pagano's  solution. 

Figure  16  compares  the  results  obtained  for  with  those  given  by  Pagano  [1978]. 
The  results  for  N=6  matched  well  with  those  of  Pagano's  except  over  the  edge 
element.  Notice  that  a  singularity  exists  at  the  edge  due  to  the  different  constitutive 
relations  at  the  free  edge.  As  the  were  interpt)lated  over  the  edge  element,  the 

edge  element  should  be  kept  relatively  small  to  better  approximate  the  singularity. 
Hence  further  spatial  refinement  was  done  by  using  15x11  and  19x11  element  meshes. 
These  results  are  shown  in  Figure  17  and  Figure  18.  These  figures  show  that  with 
mesh  refinement  near  the  edge  of  the  laminate,  a  sharp  peak,  was  observed  in  the  edge 
element.  The  magnitude  of  the  peak  increased  with  decreasing  width  of  the  edge 

element.  It  appears  that  if  the  edge  element  is  sufficiently  small  and  is  ignored,  the 

stresses  in  the  penultimate  element  can  be  considered  to  realistically  represent  the 

conditions  very  close  to  the  free  edge. 

Equation  (275)  was  used  to  compute  the  normal  stress.  The  derivatives  were 
computed  at  each  location  using  a  central  difference  method.  The  size  of  the  interval 
used  was  0.05  of  the  element  width,  reducing  to  0.02  of  the  element  width  near  the 
free  edge. 

The  results  for  cr„  for  the  different  element  meshes  are  shown  in  Figure  19. 
The  hump  noticeable  in  Pagano's  solution  was  not  observed.  It  was  observed  that  the 
edge  values  of  the  stresses  approached  Pagano’s  solution  as  the  mesh  near  the  laminate 
edge  was  refined. 
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Figure  16'.  Distribution  of  YZ-strcss  Along  0/9()  Interface  (11x11  Meshf. 
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f  igure  17:  Distribution  of  YZ-stress  Along  0/90  Interface  (15x11  Mesh). 
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Figure  18:  Distribution  of  YZ-stress  Along  0/9()  Interface  (19x11  Mcshl. 
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The  distribution  of  the  tranverse  inplane  displacement  is  shown  in  Figure  20  and 
agrees  well  with  Pagano's  result. 

6.4  NATURAL  FREQUENCIES  OF  FREE  VIBRATION  OF  A  SANDWICH 
PLATE. 

Though  the  constitutive  relations  for  transverse  shear  have  been  obtained  from  the 
equations  of  static  equilibrium,  inertial  terms  were  included  in  the  finite  element 
formulation  to  compute  the  natural  frequencies  of  free  vibration  of  a  simply  supported 
plate.  A  sandwich  plate  whose  3-D  elasticity  solution  for  the  fundamental  frequency 
was  presented  by  Srinivas  [1970]  was  considered.  The  orthotropic  laminated  plate 
consisted  of  three  layers  as  shown  in  Figure  21.  The  top  and  bottom  layers  had  the 
same  thickness  and  material  properties  while  the  thickness  and  material  properties  <rf 
the  middle  layer  were  different.  Square  geometry  with  thickness/side  ratio  of  0.1 
(h/a=h/b=0.l)  was  used  and  three  different  cases  of  material  properties  given  in  Table 
14  were  considered. 

For  a  simply  supported  laminated  plate,  the  boundary  conditions  are  [Srinivas 


1970,1973]: 

U|=H’  =  0  at  x^^Oandb  (282) 

U2  =  w  =  0  at  x^=0  and  a  (283) 

With  the  present  formulation,  these  boundary  conditions  were  restated  as 

u'^  =  <t)f  =  w  =  0  or  x^  =  0  and  b  (284) 

«y’  =  02^  =  »v  =  O  al  x^~0  and  a  (285) 


Figure  22  shows  the  finite  element  mesh  along  with  the  boundary  conditions. 


Table  14 


Lamination  Data  for  Sandwich  Plate 


Case 

tfh 

tjh 

tfh 

py ) 

1 

§2/e2 

I 

0.1 

00 

0 

0.1 

1.0 

1.0 

1 

1 

II 

0.1 

0 

00 

0.1 

1.0 

1.0 

10 

10 

III 

0.1 

0 

00 

0.1 

1.0 

1.0 

50 

50 

*  For  all  layers,  ratios  of  orthotropic  elastic  constants 
were : 

:  G22 :  G44 :  ■  G«6  =  3-802  : 0.879 : 1.9% :  1.015  : 0.608 :  1.0 


1.^4 


f  igure  22: 


I'iniie  Hlement  Mesh  and  Boundary  Conditions 


Table  15  shows  the  nondimensionalized  fundamental  natural  frequency  obtained 
using  4,  16  and  36  element  mesh  for  cases  1,  II,  III  for  which  the  ratio  of  elastic 
constants  of  the  middle  layer  to  those  of  the  outer  layers  was  1,  10  and  50, 
respectively.  It  was  seen  that  the  accuracy  for  Case  I  was  excellent,  but  the  error 

became  larger  as  the  ratio  of  the  elastic  constants  of  the  outer  layers  to  the  middle 

layer  increased.  Figure  23  illustrates  the  convergence  of  the  numerical  solutions  with 
spatial  mesh  refinement.  It  was  seen  that  convergence  after  refinement  beyond  16 
elements  was  slower  than  from  4  to  16  elements. 

i'o  examine  the  effect  of  consistent  treatment  of  transverse  shear  deformation,  the 
same  example  problem  was  solved  using  the  code  with  a  simple  shear  correction  factor 
k.=5/6,  assuming  the  constitutive  equations  of  shear  in  each  layer  to  be  uncoupled. 

Table  16  shows  the  non-dimensionalized  fundamental  natural  frequency  obtained  based 

on  this  approach  for  the  three  cases.  The  quantity  in  parentheses  is  percentage  error. 

The  results  obtained  with  k.=5/6  were  compared  with  the  result  obtained  with  the 
consistent  theory.  These  comparisons  for  cases  I,  II  and  III  are  illustrated,  respectively, 
in  Figure  24,  Figure  25  and  Figure  26.  It  was  observed  that  the  uncoupled  approach 
overpredicts  the  natural  frequency.  The  difference  increases  as  the  difference  in  stiff- 
nes.ses  of  the  outer  and  the  inner  layer  grows.  Further,  Figure  24  through  Figure  26 
show  that  the  best  accuracy  was  obtained  with  the  4x4  element  mesh,  but  the  solution 
did  not  show  monotonic  convergence. 

To  further  study  the  effect  of  the  proposed  consistent  shear  constitutive  relations, 
in  addition  to  comparing  the  fundamental  frequency  the  higher  frequencies  obtained 
using  the  proposed  theory  and  a  shear  correction  factor  of  5/6  were  also  compared. 
Figure  27  to  Figure  29  compare  the  frequencies  of  a  537  degree-of-freedom  system 
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Table  15 


Non-Dimensionalized  Fundamental  Frequency  by  FEM  based  on  the  Con¬ 
sistent  Shear  Deformable  Theory 


Mesh 

Case  1 

Case  II 

Case  III 

4 

0.094697 

0.19682 

0.31097 

(2.4%) 

(2.9%) 

(3.8%) 

16 

0.092952 

0.19472 

0.30919 

(0.5%) 

(1.8%) 

(3.2%) 

36 

0.092900 

0.19463 

0.30911 

(0.4%) 

(1.7%) 

(3.1%) 

Exact  3-D 

0.09248 

0.19132 

0.29954 

(Srinivas) 

•  X  = 


where  (O  is  natural  frequency. 
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PERCENT  ERROR 


Table  16 


Non-Dimensionalized  Fundamental  Frequency  by  FEM  using  shear  correc¬ 
tion  factor  k=5/6 


Mesh 

Case  I 

Case  II 

Case  III 

4 

0.09498 

0.19710 

0.31489 

(2.7%) 

(3.0%) 

(5.1%) 

16 

0.09317 

0.19584 

0.31305 

(0.7%) 

(2.3%) 

(4.5%) 

36 

0.09341 

0.19805 

0.32232 

(0.7%) 

(3.5%) 

(7.6%) 

Exact 

0.09248 

0.19132 

0.29954 

(Srinivas) 
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PROPOSED  THEORY  0 
WITH  K=5/6  A 


o 

o 


ligure  24:  Lffect  of  (Consistent  Shear  (Correction  for  Case  1 
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32 


Figure  25:  Effect  of  Consistent  Shear  Correction  for  Case  11 
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PERCENT  ERROR 

-10.00  -5.00  0.00  5.00  10.00 


PROPOSED  THEORY  0 
WITH  K=5/6  A 


2 


U  8  16 

NO.  OF  ELEMENTS 


32 


Figure  26: 


{Effect  of  Consistent  Shear  Correction  for  (Tase  III 


64 
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(4x4)  mesh  for  the  three  cases  of  the  sandwich  plate  studied.  The  values  of  the  lower 
frequencies  computed  using  the  consistent  shear  approach  differ  only  slightly  from 
those  by  the  discrete  laminate  using  a  shear  correction  factor  of  5/6.  However,  the 
higher  frequencies  are  distinctly  different. 
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FREQUENCY 


Section  VII 


DISCUSSION 

I^ue  to  the  inherent  complexity  associated  with  material  anisotropy  and 
tnhomogeneity  of  composites,  a  laminated  plate  often  shows  quite  different  mechanical 
characteristics  from  the  homogeneous  isotropic  counterpart.  Therefore,  it  is  essential 
that  any  simplified  theory  satisfy  equilibrium,  kinematic  and  constitutive  relations  as 
closely  as  possible  to  ensure  reliable  results. 

The  investigatitm  reported  here,  aimed  at  development  of  stress  and  deformation 
analysis  of  laminated  composites  resulted  in  the  following  accomplishments: 

a.  Development  of  a  systematic  and  general  theory  to  consistently  incor¬ 
porate  the  transverse  shear-deformation  effect  in  composite  laminates. 

b.  Derivation  of  variational  principles  governing  the  refined  theory  to 
provide  a  basis  for  development  of  efficient  and  reliable  Ritz  type  as 
well  as  finite  element  approximation  procedures. 

c.  Implementation  of  the  theory  into  a  finite  element  computer  program 
including  code  verification. 

d.  Application  of  the  theory  to  free-edge  delamination  specimens. 

A  discrete  laminated  plate  theory  based  on  the  assumption  of  a  layerwise  linear 
variation  of  in-plane  displacements  has  been  further  refined  by  incorporating  the  effect 
of  transverse  shear  deformation  in  a  consistent  manner,  viz.  allowing  for  the  coupling 
of  shear  deformation  effects  between  layers.  This  development  was  facilitated  by  a 
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mixed  variational  principle  of  linear  elasticity  derived  using  Sandhu's  generalized 
prcKedure  for  variational  formulation  of  linear  coupled  boundary  value  problems.  This 
mixed  variational  principle  is  more  useful  for  the  application  to  a  general  anisotropic 
material  than  Reissner's  [1984]  approach.  A  parabolic  distribution  of  shear  stresses  over 
the  thickness  of  each  layer  was  assumed.  Continuity  of  stresses  and  displacements  in 
the  layer  interfaces  were  allowed  for.  Distinctive  features  of  the  resulting  constitutive 
relations  for  transverse  shear  are: 

1.  Shear  force  resultants  for  each  layer  are  a  linear  combination  of  the  shear 

strains  of  all  the  layers.  Directional  coupling  of  the  constitutive  relations  dis¬ 
appears  for  orthotropically  constructed  laminates. 

2.  Coefficients  in  a  linear  combination  of  shear  strains  are  determined  by  parame¬ 
ters  related  to  lamination  schemes  such  as  material  properties,  thickness  of  lay¬ 

ers  and  stacking  sequence  of  layers. 

3.  The  shear  constitutive  relations  also  depend  upon  tangential  stresses  specified 

on  the  laminate  surfaces. 

The  fact  that  the  shear  force  over  a  layer  is  coupled  with  the  shear  strains  of  other 
layers  in  a  linear  fashion  seems  to  be  a  striking  result  because  such  coupling  has  not 
been  anticipated  in  earlier  work.  However,  this  result  can  be  attributed  to  the  consid¬ 
eration  of  continuity  of  shear  stresses  in  the  interface  of  layers  which  has  not  been 
taken  into  account  in  previous  studies.  Coupling  of  shear  constitutive  relations  of  all 
layers  could  result  in  a  better  understanding  of  how  a  laminate  composed  of  many 
layers  might  react  to  applied  forces. 

The  nature  and  extent  of  the  shear  coupling  was  studied  by  looking  at  the 
constitutive  matrix  for  a  12-layer  graphite  epoxy  laminate.  It  was  observed  that  the 
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shearing  force  in  a  layer  was  not  significantly  influenced  by  the  strains  in  other 
layers.  Also,  the  effect  was  local,  i.e,  contribution  from  other  layers  decreased  sharply 
with  increasing  distance  between  layers. 

Using  the  refined  laminate  theory,  a  systematic  development  of  variational 

principles  for  static  as  well  as  dynamic  analysis  of  laminated  composite  plates  was 
carried  out.  Direct  as  well  as  complementary  formulations  were  developed.  The 
complementary  self-adjoint  form  of  the  field  equations  obtained  is  an  advancement  over 
the  one  presented  by  .AI-Ghothani  [1986],  insomuch  as  the  present  work  contains 
coupling  terms  of  the  transverse  shear  con.stitutive  equations  between  layers. 
Nonhomogeneous  boundary  conditions  and  internal  jump  discontinuities  have  been 
explicitly  included  in  the  general  variational  principles.  Allowance  for  jump 
discontinuity  terms  in  the  variational  formulation  is  meaningful  in  the  context  of 
direct  approximation  in  finite  element  spaces  since  the  space  of  approximants  may  not 
be  sufficiently  smooth.  Extensions  of  the  general  variational  principles  through 

elimination  of  certain  field  operators  and  specializations  by  restricting  some  of  the  field 
equations  and/or  boundary  conditions  to  be  identically  satisfied  have  been  proposed. 

Based  on  a  special  variational  principle,  a  finite  element  procedure  which  uses 

u  and  w  as  the  field  variables  has  been  formulated  and  a  finite  element  code 

has  been  developed.  The  computer  program  was  used  to  study  the  effect  of  the 
constitutive  coupling  by  solving  the  static  problem  of  a  laminated  coupon  under  axial 
extension.  The  results  were  compared  with  the  solution  provided  by  Pagano  [1978]. 
The  displacement  solution  was  seen  to  agree  well.  In  addition  to  refinement  of  the 

Finite  Element  mesh,  increased  accuracy  was  seen  to  result  from  the  division  of  each 
layer  into  sublayers.  The  results  for  inplane  stresses  compared  well  though  those  for 
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transverse  shear  and  normal  stresses  did  not.  It  was  observed  that  refinement  of  the 
finite  element  mesh  along  the  length  of  the  free-edge  delamination  specimen  did  not 
significantly  contribute  to  improvement  in  accuracy.  However,  refinement  in  the  lateral 
direction  (y-direction  or  the  axis)  near  the  free-edge  gave  considerable  improvement. 
Also,  refinement  of  layers  into  sublayers  improved  the  results.  This  points  to  the 
desirability  of  using  a  higher  order  variation  over  the  x,  co-ordinate. 

Though  the  constitutive  relations  have  been  derived  from  the  equations  of  static 
equilibrium,  the  problem  of  free  vibration  of  a  sandwich  plate  was  also  studied  by 
including  inertial  terms  in  the  variation  formulation.  The  fundamental  frequency  was 
compared  with  that  obtained  from  an  elastodynamic  solution.  It  was  observed  that 

1.  In  certain  cases,  a  consistent  shear  correction  improves  the  accuracy  considera¬ 
bly  in  predicting  natural  frequencies  over  the  shear  correction  by  a  simple 
factor  which  has  been  widely  used  in  previous  theories. 

2.  Improvement  of  accuracy  depends  upon  the  material  properties  of  each  layer. 
With  increase  in  the  difference  of  material  properties  in  the  individual  layers, 
the  significance  of  a  consistent  shear  correction  was  more  pronounced. 

The  higher  frequencies  were  also  studied  and  it  was  observed  that  the  higher 
HKxles  differed  sharply  f rom  those  obtained  using  a  single  shear  correction  factor. 

These  limited  numerical  tests  illustrate  the  validity  of  a  consistent  shear  correction 
procedure.  Considering  numerous  design  possibilities  of  composite  laminates,  use  of  a 
procedure  to  treat  transverse  shear  deformation  in  a  consistent  manner  rather  than 
shear  correction  on  an  ad-hoc  basis  seems  to  be  desirable. 
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Appendix  A 

VARIATIONAL  FORMULATION  OF  LINEAR  COUPLED 

PROBLEMS 


The  procedure  for  obtaining  variational  principles  governing  linear  coupled 
boundary  value  problem  is  summarized  following  Sandhu  [1970,1971,1975,1976].  I'lie 
procedure  can  be  considered  as  an  extension  of  .MikhJin's  [1965]  basic  variational 
theorem  to  coupled  linear  boundary  value  problems,  including  nonhomogeneous  boundary 
conditions  and  internal  discontinuities  which  may  exist  in  the  physical  problem  or  arise 
in  connection  with  numerical  approximation  procedures. 

A.1  MAT/ll-MATICAL  PRELIMINARIES 

A.  1.1  Boundary  Value  Problem 

Consider  the  boundary  value  problem 

.4(uJ  =  /  on  R  (A.l) 

Ciu)  =  g  on  aR  (A.2) 

here  A  and  C  are  the  linear,  bounded  opierators,  u  is  the  field  variable,  R  is  an  open 
connected  region  in  an  Euclidean  space  and  a^  is  i^s  boundary.  Let  and  be 

linear  vector  spaces  defined  on  the  regions  indicated  by  subscripts  such  that 
/€V^  and  g€V^^ 
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I'hen.  the  operators  can  be  regarded  as  transformations  defined  over  sets  such 

that 

A:W^^V^  (A.3) 

F'or  A,  C  differential  operators,  are,  in  general,  dense  subsets  in  V’^, 

respectively. 

A. 1.2  Bilinear  Mapping 

I.et  and  S  be  linear  vector  spaces.  A  bilinear  mapping  B(w  ,v);  \  x\'  — *  <S’  assigns 
an  element  in  S  to  an  ordered  pair  of  elements  w,i>6\'  while  preserving  linearity,  lor 
convenience,  we  shall  use  the  notation 

B(w,v)  =  <>v,v>  (A.5) 

Bilinear  mapping  B  is  said  to  be  nondegeneraie  if 

<w,v>  =0  for  all  w  if  and  oniy  if  v  =  0  (.^.6) 

and  symmetric  if 

<w,  v>  =  <v,w>  for  all  v,w  €V^  (A.7) 

A.1.3  Self-Adjoint  Operator 

I.et  — \'  be  an  operator  on  the  linear  vector  space  V.  Operator  is 

said  to  be  adpint  of  A  with  respect  to  a  bilinear  mapping  <  ,  >^:  VxV'-»,S,  a  linear 
vector  space,  if 

<w,  Av>  =  <v,  a’w>ji  +  D^jfv.w)  (A.8) 

for  all  v,w  €  V.  We  assume  here  that  the  domain  of  A,  A\  a  dense  subset  in  V  can  be 
extended  to  V.  Here,  the  subscript  R  of  bilinear  mapping  indicates  that  V  is  defined 
over  spatial  region  R  and  Dj,/v,w)  represents  quantities  as.sociated  with  the  boundarv 
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QR  of  R.  If  =  A\  a  is  said  to  be  self-adjoint.  In  particular,  an  operator  A  is 
said  to  be  symmetric  with  respect  to  bilinear  mapping  if 

<w ,  Av> =  <v,  Aw> g  (A.9) 

A. 1.4  Gateaux  Differential 

Consider  a  continuous  function  fl  :  V -*  5.  Gateaux  differential  of  fl  is  defined  as 

A  fl(u)  =  lim [fl(ii -I- Av)  —  n(u)]  (A.IO) 

'■  K^O  A 

provided  the  limit  exists.  Here,  v  is  referred  to  as  the  'path'  and  A  is  a  scalar.  We 
note  that  for  u ,  v€\  ,  u  +  Xv  €.  V.  If  the  Gateaux  differential  exists  at  every  point  in  a 
neighborhood  of  u  = 

A  ft(u)  =  -^fl(u-fAu)i 

V  \-u 

A.2  BASIC  VARIATIONAL  PRINCIPLE 

For  the  boundary  value  problem  (A-l)  with  homogeneous  boundary  condition, 
Mikhlin  [1965]  used  the  inner  product  as  the  nondegenerate  symmetric  bilinear  mapping 
on  the  linear  vector  space  of  square  integrable  functions  and  showed  that  the  unique 
solution  minimizes  the  functional 

fl(u)  =  <  Au ,  u>  ^  —  2<u ,  f>  (A. 11) 

if  the  linear  operator  A  is  positive  definite  and  self-adjoint.  Conversely,  which 

minimizes  the  functional  (A-11)  is  the  solution  of  the  problem  (A-l). 

Taking  Gateaux  differential  of  the  function  (A-ll), 

A  fl(u)  =  lim  -^  [<  Aiu  +  Av>,u  +  Av>  — 2<u  +  Av>  <  Au ,  u>  +2<u,/>] 

'  \-o  A 


=  <Au,v>  +  <Av,u>  — 2<v,/> 


=  2  <  V ,  Au~  f> 


(A. 12) 


The  Gateaux  differential  evidently  vanishes  at  the  solution  u  =  such  that 
Au^  —f  =  0.  Also,  since  the  bilinear  mapping  is  nondegenerate,  vanishing  of  the 
Gateaux  differential  for  all  v  implies  Au^  —  f  =  0.  To  prove  the  minimization  property, 
the  bilinear  mapping  has  to  be  into  the  real  line  and  the  operator  must  be  positive. 
However,  in  general,  it  is  only  necessary  to  use  vanishing  of  Gateaux  differential  as 
equivalent  to  (A-1)  being  satisfied. 

AJ  VARIATIONAI,  lORMVLAT/ON  OF  TUK  COVPLED  PROBLEM 

Through  generalization  of  Mikhlin's  theorem,  Sandhu  [1970,1971,1975,1976] 
constructed  a  framework  to  handle  the  inverse  problem  of  variational  calculus  for  the 
linear  coupled  problem  with  multiple  field  variables.  To  include  the  nonhomogeneous 
boundary  conditions  and  internal  discontinuity  conditions  explicitly  in  the  formulation, 
the  concept  of  consistent  boundary  operators  was  introduced. 

Consider  the  coupled  boundary  value  problem  with  multiple  field  variables 

A(u)  =  f  on  R  (A. 13) 

C(u)  =  g  on  QR  (.A. 14) 

in  which  A ,  C  and  u  are.  respectively,  the  field  operator  matrix,  the  boundary  operator 
matrix  and  the  vector  of  field  variables,  and  f,g  are  the  vectors  of  known  forcing 

functions.  Explicitly, 

n 

La  u  =  f  on  R  (A. 15) 

r  > 

=  on  a/2,.  i=  1,  2 . n  (A. 16) 
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in  wliith  QK  denote  segments  of  QK  such  that 


e/e  =  U  (A.17) 

1=1 

and  n  is  the  number  of  independent  field  variables.  Operators  are  regarded  as  the 
transformations 

.4  :  M  -  P  (A. 18) 

v  here 

•? 

u  e  V\-  =  n  (A. 19) 

I 


n 

y=t 

Thus,  the  range  of  A,^  constitutes  a  product  space 

V  =  \\x\\x . V'  (A.21) 

12  n 

I.et  V  be  a  vector  space  defined  as  the  direct  sum 

V  =  V\+V,  + . +  V  (A.22) 

12  n 

and  an  element  u€V  be  the  ordered  set 

u  =  }  (A. 23) 

1  ^  r. 


such  that  «, €\',.  Then,  a  bilinear  mapping  on  V',  may  be  defined  as 

<u ,  v>  .,  =  <t/, ,  V,  >  „  + . +  <u  ,  V  >  „  (A.24) 

A’  I  I/?  n  n  K 

The  set  of  operators  A,^  is  said  to  be  self-adjoint  with  respect  to  this  bilinear  mapping 
if  (Sandhu,  1976] 
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\'>here  ,u,)  denote  quantities  associated  with  ^R.  If  the  set  of  operators  A  is 

self  adpint.  as  a  generalization  of  Mikhlin's  theorem,  the  function  governing  the 
problem  (A-13)  and  (A-14)  was  defined  as 

n  n 

n  =  +  <u,,Cu-2g>^;i  (A.26) 

.■1  r\ 

Tor  vanishing  of  the  Gateaux  differential  of  this  function  to  imply  (A-13)  and  (A-14). 
the  Ixtundary  operat(irs  must  be  consistent  with  the  field  operators  A,.  Sandhu 
[1076]  stated  the  consistency  condition  of  the  Ixuindarv  operators  as 

n  r/ 

/  1  /I 

In  other  words,  for  (A'26)  to  be  the  governing  function  in  variational  formulation  of 
the  problem  given  by  (A-13)  and  (A-14),  the  boundary  operators  must  satisfy  (A-27). 
Sandhu  [1975]  showed  that  appropriate  boundary  terms  should  be  included  in  the 
formulation  even  if  the  boundary  conditions  are  homogeneous. 

In  applying  the  finite  element  method  to  obtain  an  approximate  solution,  the 
approximation  space  may  have  limited  smoothness,  e.g.,  when  nonconforming  elements 


are  used,  Prager  [1968] 

was 

the  first 

to  discuss  this 

aspect 

in  connection  with 

variational  formulation. 

1o 

allow  for 

this,  Sandhu 

[1975] 

wrote 

the  internal 

discontinuity  conditions  in 

the 

form 

(C  u  y  =  g' 

on 

(A.28) 

where  a  prime  denotes 

the 

internal  jump  discontinuity 

along 

internal 

surfaces  Q/J, 

embedded  in  the  region  R.  Since  (A-28)  has  the  same  form  as  the  boundary 
conditions,  it  is  apparent  that  this  condition  can  be  included  in  the  governing  function 
in  the  same  way.  The  governing  function  allowing  for  (A-28)  is 
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(A. 29) 


1 


1 


/■I  /-I 


This  is  the  general  form  of  the  governing  function  in  the  variational  formulation  of 
linear  coupled  boundary  value  problem  with  multiple  field  variables.  The  essential 
step  in  setting  up  a  variational  formulation  of  a  boundary  value  problem  is  to  write 
the  field  eouations  in  a  form  that  the  matrix  of  field  operators  is  self-adpint  with 
respect  to  eertair  bilinear  mapping  and  the  boundary  conditions  are  consistent  with  the 
held  o]-)erators.  1  he  procedure  is  also  applicable  to  initial-boundary  value  problems  using 
appropriate  bilinear  mappings. 


0 
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Appendix  B 

EVALUATION  OF  ELEMENT  MATRICES 

B.l  INTKRPOLATION  FUNCTIONS  OF  THE  ‘HETEROSIS  ELEMENT 

The  computer  program  developed  incorporated  the  'Heterosis'  plate  Irending  dement 
(Hughes  1978]  along  with  reduced/selective  integration  technique.  The  element  matrices 
can  be  formed  following  the  usual  procedure  of  isoparametric  element  formulation. 
However,  the  'Heterosis'  element  differs  from  other  isoparametric  elements  in  using 
different  interpolation  scheme  for  lateral  displacement  on  the  one  hand  and  inplane 
displacement  and  rotations  of  the  cross-section  on  the  other.  In-plane  displacements 
and  rotations  of  cross-section  are  approximated  by  qaudratic  functions  for  8-node 
isoparametric  element  while  the  lateral  displacement  w  is  approximated  by  9-node 
l.agrange  interpolation  functions.  Consequently,  the  number  of  degrees  of  freedom  at 
the  center  node  is  less  than  that  at  other  nodes  by  one.  Interpolation  scheme  of  the 
Heterosis'  element  is  showm  in  Figure  B.I.  Using  such  interpolation  scheme,  Hughes 
[1978]  was  able  to  avoid  spurious  zero  energy  mode  of  stiffness  matrix  which  can  be 
caused  by  use  of  reduced  integration. 

Interpolations  functions  of  8-node  isoparametric  element  and  9-node  Lagrange 
element  in  terms  of  natural  coordinates  (s,t)  and  their  derivatives  with  respect  to  s 
and  t  are  given  below. 
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(1-sXl-rX-l-s-i) 

(l-tX2s-t-t) 

(1-5x21-1-5) 

(i-hsXl-rX-l+5-t) 

(i-fX25-r) 

(1-1^5x21-5) 

(l+sXl-t-tX-l-t-s+t) 

(l+rX25-(-i) 

(1-1-5x21-1-5) 

1 

(l-sXl+tX-l-s+r) 

6N  _  1 

(l-t-rX25-/) 

aN  _  1 

(1-5x21-5) 

4 

2(l-j’Xl-r) 

QJ  4 

-4  5(1-0 

4 

-2(1-5") 

2(l+sXl-t') 

2(l-t') 

-4  Kl-i-s) 

2(l-s’Xl+r) 

-4  j(l+r) 

2(1-5") 

2(l-5Xl-t^) 

-2(l-r') 

-4 1(1-5) 

St  (1  -5)(l-l) 

i(25-1KM) 

5(21  1)(5  1) 

51  (l-isKM) 

1(25-i-1)(i-1) 

5(2i-1  )(5-1-1  ) 

5/  (  1 -1-5  )  (  1 -1-/ ) 

l(25-l-l)(l-hl) 

5(2/-i1)(5+1) 

St  (5-I  )(l-l  1  ) 

i(25-1)(i-i^1) 

5(2/+1)(5-1) 

2i(1-5")(i-1) 

it  =  1 

451  (1-1) 

it  =  1 

2(2l-hl)(l-5") 

25(5-1-1)(1-i") 

as  4 

2(25+1)(1-i") 

ai  4 

-4  51  (5-1-1) 

2i(1-5")(i+1) 

-451(1-1-1) 

2(2l-il)(l-5") 

25(5-1)(m") 

2(25-1)(1-i") 

451(1-5) 

4(1-5")(1-i") 

85(1"-!) 

81  (5"-!) 

Mere,  N  and  L  denote  interpolation  functions  for  8-node  isoparametric  and  9-node 


Lagrange  element:,  respectively. 
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R.2  EVALUATION  OF  STIFFNESS  AND  MASS  MATRICES 

Since  the  field  variables  are  interpolated  over  an  element  in  natural  coordinates 
(s,t),  it  is  necessraj"  to  set  up  the  relation  of  the  global  coordinates  and  natural  (local) 
coordinates  for  evaluation  of  the  element  matrices  defined  in  Section  VI.  We  consider 
a  mapping  of  global  coordinate  system  (jt,  to  local  coordinate  system  (s,t).  We 
assume  that  this  mapping  is  one-to-one  and  onto.  By  chain  rule,  the  derivative  in 
each  coordinate  system  is  related  by 


164 


(B.3) 


e 

9 

9 

9 

05 

=  J 

9^ 

rtf 

dx 

=  J  ‘ 

9^ 

9 

9 

Ur 

—  J 

9 

d* 

9y 

9y 

9t 

where  Jacobian  matrix  J  and  its  inverse  is  defined  as 


J  = 


9z 

iy  -ill 

9^  9^ 

and  r*  =  J- 

9* 

Qs 

Bx  By 

IJI 

-i£  i£ 

Bt  Bf 

9* 

Br  , 

(B.4) 


Here,  IJI  is  the  determinant  of  Jacobian  matrix.  Using  (B-3)  and  (B'4),  one  can  obtain 
the  expressions  of  the  matrices  T^,  and  defined  in  (214)-(216)  in  natural 
coordinates,  hollowing  the  concept  of  isoparametric  formulation,  global  coordinates  are 
interpolated  over  an  element  as 


x  =  'i'  X 


(B.5) 


where  'If  is  the  vector  of  interpolation  functions  used  for  field  variable,  x  is  the 
vector  of  global  coordinate  values  at  nodal  points.  Since  in  the  Heterosis'  element 
different  interpolations  functions  N  and  L  are  used  for  0“^  and  w,  ^  in  (B-5)  must  be 
L  for  evaluating  T„  and  while  be  N  for  Tj. 

Using  (B-IMB-5)  and  the  interpolation  functions  defined  in  (21l)-(216), 


H 


L  0 
0  L 


(B.6) 


H  =  N 

b 


(B.7) 


= 


1 


y.  0 
0 

r/  / 


w  here 


(R8) 


(R9) 
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IJJ  =  -/Pi! 


ij^i  = 

R  =  N,N.f-N.,N.f 


Here,  a  subscripted  comma  denotes  partial  differentiation  with  respect  to  the  variables 
following  the  comma. 

In  element  matrices  given  in  (218)-(233),  the  integrands  are  functions  of  natural 
coordinates  (s-t).  Therefore,  the  surface  integration  extends  over  the  natural  coordinate 
surface.  Since,  in  general. 


dR  =  \J\d5dt 

integration  in  each  coordinate  system  is  related  by 


1  1 


J*  j Hxjdxdy  =  J  jFisjdmsdt 


-1  -1 


Using  Gaussian  quadrature 


(B.10) 


(B.11) 


f  F(x,y)dR  = 

'y  ;-l  y-I 

where  m  is  the  number  of  Gaussian  quadrature  points  and  W,^  are  weighting  values. 
Here,  it  should  be  mentioned  that  the  in  the  'Heterosis'  element  numerical  integration 
was  performed  by  selective/reduced  integration  technique,  viz.,  two-point  Gaussian 
quadrature  for  evaluation  of  transverse  strain  energy  term  while  three-point  quadrature 
is  used  for  other  quantities.  Therefore,  (231M233)  were  evaluated  by  two-point 
quadrature  and  the  remaining  quantities  were  evaluated  by  three-point  quadrature. 
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